{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "DontBlameTheELBO.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MhWvn9RtrOZQ",
        "colab_type": "text"
      },
      "source": [
        "Copyright 2019 Google LLC.\n",
        "\n",
        "Licensed under the Apache License, Version 2.0 (the \"License\");"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OcoAHDwJOrOQ",
        "colab_type": "text"
      },
      "source": [
        "# The Linear VAE\n",
        "\n",
        "### The model:\n",
        "\n",
        "$\n",
        "\\newcommand{\\bz}{\\mathbf{z}}\n",
        "\\newcommand{\\bx}{\\mathbf{x}}\n",
        "\\newcommand{\\bu}{\\mathbf{u}}\n",
        "\\newcommand{\\bI}{\\mathbf{I}}\n",
        "\\newcommand{\\bW}{\\mathbf{W}}\n",
        "\\newcommand{\\bD}{\\mathbf{D}}\n",
        "\\newcommand{\\bC}{\\mathbf{C}}\n",
        "\\newcommand{\\bS}{\\mathbf{S}}\n",
        "\\newcommand{\\bV}{\\mathbf{V}}\n",
        "\\newcommand{\\bL}{\\mathbf{L}}\n",
        "\\newcommand{\\bU}{\\mathbf{U}}\n",
        "\\newcommand{\\bR}{\\mathbf{R}}\n",
        "\\newcommand{\\bLambda}{\\boldsymbol{\\Lambda}}\n",
        "\\newcommand{\\bmu}{\\boldsymbol{\\mu}}\n",
        "\\newcommand{\\cN}{\\mathcal{N}}\n",
        "\\newcommand{\\reals}{\\mathbb{R}}\n",
        "\\newcommand{\\expect}{\\mathbb{E}}\n",
        "$\n",
        "Suppose $\\bz \\in \\reals^q$ generates $\\bx \\in \\reals^d$. We use a standard Gaussian prior:\n",
        "\n",
        "\n",
        "$$p(\\bz) = \\cN(\\mathbf{0},\\bI)$$\n",
        "\n",
        "And define a linear generative model from the latent space with a Gaussian observation model.\n",
        "\n",
        "$$p(\\bx|\\bz) = \\cN( \\bW \\bz + \\bmu, \\sigma^2 \\bI)$$\n",
        "\n",
        "This model is known as [Probabilistic PCA (pPCA)](http://www.robots.ox.ac.uk/~cvrg/hilary2006/ppca.pdf). A linear VAE can be used to maximize the likelihood of pPCA.\n",
        "\n",
        "### Optimal pPCA MLE parameters\n",
        "\n",
        "The maximum likelihood estimate of $\\bmu$ is the mean of the data. We can compute $\\bW_{MLE}$ and $\\sigma^2_{MLE}$ as follows:\n",
        "\n",
        "$$\\bW_{MLE} = \\bU_{q}(\\bLambda_q - \\sigma^2 \\bI)^{1/2}\\bR,\\\\\n",
        "\\sigma^2_{MLE} = \\frac{1}{n - q} \\sum_{j=q+1}^{n}\\lambda_j$$\n",
        "\n",
        "Here $\\bU_q$ corresponds to the first $q$ principal components of the data with the corresponding eigenvalues $\\lambda_1, \\ldots, \\lambda_q$ stored in the $q\\times q$ matrix $\\bLambda_q$. The matrix $\\bR$ is an arbitrary rotation matrix which accounts for weak identifiability in the model. The MLE for $\\sigma^2$ corresponds to the average variance lost in the projection.\n",
        "\n",
        "## Learning pPCA with a Linear VAE\n",
        "\n",
        "We will aim to learn a variational distribution $q(\\bz|\\bx)$ to approximate the posterior $p(\\bz|\\bx)$. The variational bound is characterized as follows:\n",
        "\n",
        "$$\n",
        "\\log p(\\bx) = D_{KL}(q(\\bz|\\bx)||p(\\bz|\\bx)) + \\expect_{q(\\bz|\\bx)}[−\\log q(\\bz|\\bx) + \\log p(\\bx,\\bz)], \\\\\n",
        "\\Rightarrow \\log p(\\bx)  \\geq \\expect_{q(\\bz|\\bx)}[−\\log q(\\bz|\\bx) + \\log p(\\bx,\\bz)], \\\\\n",
        "= −D_{KL}(q(\\bz|\\bx)||p(\\bz)) + \\expect_{q(\\bz|\\bx)}[\\log p(\\bx|\\bz)] \\qquad (:= ELBO)\n",
        "$$\n",
        "\n",
        "Let us fix an amortized, linear Gaussian variational distribution:\n",
        "\n",
        "$$q(z|x) = \\cN(\\bV(\\bx - \\bmu), \\text{diag}(\\sigma^2_1, \\ldots, \\sigma^2_q))$$\n",
        "\n",
        "We use a diagonal covariance matrix (denoted $\\bD$ later) which is shared across all data points. We show in our paper that this is sufficient to represent the true posterior distribution. Thus, the linear VAE corresponds exactly to fitting a pPCA model with variational inference.\n",
        "\n",
        "### Exact ELBO computation \n",
        "\n",
        "The ELBO is made up of two components, the prior KL divergence (i), $KL(q(\\bz|\\bx)||p(\\bz))$ and the reconstruction loss (ii), $\\expect_{q(\\bz|\\bx)}[\\log p(\\bx|\\bz)]$.\n",
        "\n",
        "Because everything is Gaussian, we can compute these terms in closed-form.\n",
        "\n",
        "(i): $KL(q(\\bz|\\bx)||p(\\bz)) = 0.5(tr(\\bD) -\\log\\det(\\bD) + (\\bx-\\bmu)\\bV^T\\bV(\\bx - \\bmu) - q)$\n",
        "\n",
        "(ii): $\\expect_{q(\\bz|\\bx)}\\left[ \\log p(\\bx|\\bz) \\right] = \\frac{1}{2\\sigma^2} [ -tr(\\bW\\bD\\bW^T) - (\\bx - \\bmu)^T\\bV^T\\bW^T\\bW\\bV(\\bx-\\bmu) + 2(\\bx - \\bmu)^T\\bW\\bV(\\bx-\\bmu) - (\\bx - \\bmu)^T(\\bx - \\bmu)] - \\frac{d}{2}\\log 2\\pi\\sigma^2$\n",
        "\n",
        "Thus, we can train the linear VAE using exact gradients to validate the theoretical results presented in our paper."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "c-RergnbOpqp",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "%tensorflow_version 1.x\n",
        "import tensorflow as tf\n",
        "from tensorflow.data import Dataset\n",
        "import tensorflow_datasets as tfds\n",
        "\n",
        "import numpy as np\n",
        "import matplotlib as mpl\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "_LOG_2PI = np.log(2 * np.pi)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XL6pGaQ5O-HQ",
        "colab_type": "text"
      },
      "source": [
        "## Data loading utilities"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yPg6mvZLOwLt",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def load_mnist_data(take):\n",
        "  '''Load a subset of the MNIST dataset, repeated.\n",
        "\n",
        "  Args:\n",
        "    take: The amount of training data to take.\n",
        "  \n",
        "  Returns:\n",
        "    A shuffled subset of the MNIST training dataset.\n",
        "  '''\n",
        "  mnist_builder = tfds.builder(\"mnist\")\n",
        "  mnist_info = mnist_builder.info\n",
        "  mnist_builder.download_and_prepare()\n",
        "  datasets = mnist_builder.as_dataset()\n",
        "  train_dataset = datasets['train']\n",
        "  train_dataset = train_dataset.take(take).map(\n",
        "      lambda x : tf.reshape(tf.cast(x['image'], tf.float32) / 255.0, [-1]))\n",
        "  return train_dataset\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GkTPHP-DPDB9",
        "colab_type": "text"
      },
      "source": [
        "## Linear VAE"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "znbhLJqhPBFZ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class LinearVAEConfig(object):\n",
        "  \"\"\"Config class for the Linear VAE.\n",
        "\n",
        "  Args:\n",
        "    input_dim: Number of input dimensions (default: 784)\n",
        "    hidden_dim: Number of hidden dimensions (default: 20)\n",
        "    trainable_sigma: Whether observation noise should be learned (default: True)\n",
        "    sigma_init: Initial value of sigma^2 (default: None)\n",
        "    stochastic_ELBO: Use stochastic ELBO estimation (default: False)\n",
        "  \"\"\"\n",
        "\n",
        "  def __init__(self,\n",
        "               input_dim=784,\n",
        "               hidden_dim=20,\n",
        "               trainable_sigma=True,\n",
        "               sigma_init=None,\n",
        "               stochastic_ELBO=False):\n",
        "    self.input_dim = input_dim\n",
        "    self.hidden_dim = hidden_dim\n",
        "    self.trainable_sigma = trainable_sigma\n",
        "    self.sigma_init = sigma_init\n",
        "    self.stochastic_ELBO = stochastic_ELBO\n",
        "\n",
        "class LinearVAE(object):\n",
        "  \"\"\"Linear VAE\n",
        "  \n",
        "  Allows for the reconstruction error to be computed analytically or via\n",
        "  single-sample stochastic estimation.\n",
        "  \"\"\"\n",
        "\n",
        "  def __init__(self, input_tensor, config=LinearVAEConfig()):\n",
        "    self._config = config\n",
        "    self.input_tensor = input_tensor\n",
        "\n",
        "    self.z_logvar = tf.get_variable(\n",
        "        'logvars', shape=(config.hidden_dim), trainable=True)\n",
        "    self.expanded_z_logvar = tf.expand_dims(self.z_logvar, 0)\n",
        "\n",
        "    self.enc_weight = tf.get_variable(\n",
        "        'enc_weight',\n",
        "        shape=(config.input_dim, config.hidden_dim),\n",
        "        trainable=True)\n",
        "    self.mean = tf.get_variable(\n",
        "        'mean', shape=(config.input_dim), trainable=True)\n",
        "    self.dec_weight = tf.get_variable(\n",
        "        'dec_weight',\n",
        "        shape=(config.hidden_dim, config.input_dim),\n",
        "        trainable=True)\n",
        "    if self._config.sigma_init is not None:\n",
        "      sigma_init = tf.constant_initializer(np.log(self._config.sigma_init))\n",
        "    else:\n",
        "      sigma_init = None\n",
        "    self.log_sigma_sq = tf.get_variable(\n",
        "        'log_sigma',\n",
        "        shape=(),\n",
        "        initializer=sigma_init,\n",
        "        trainable=self._config.trainable_sigma)\n",
        "    self.sigma_sq = tf.exp(self.log_sigma_sq)\n",
        "\n",
        "    self.z_mean = tf.matmul(self.input_tensor - self.mean, self.enc_weight)\n",
        "    if self._config.stochastic_ELBO:\n",
        "      eps = tf.random_normal(\n",
        "        (tf.shape(self.input_tensor)[0], self._config.hidden_dim),\n",
        "        mean=0.0,\n",
        "        stddev=1.0,\n",
        "        dtype=tf.float32)\n",
        "      self.z_sample = self.z_mean + eps * tf.exp(0.5 * self.expanded_z_logvar)\n",
        "\n",
        "      with tf.variable_scope('decoder') as dec_scope:\n",
        "        self.reconstruction = tf.matmul(self.z_sample, self.dec_weight)\n",
        "\n",
        "    ## Loss\n",
        "    with tf.name_scope('losses'):\n",
        "      # Reconstruction loss under recognition model\n",
        "      if self._config.stochastic_ELBO:\n",
        "        self.recon_loss = self.stochastic_recon(self.input_tensor,\n",
        "                                                self.reconstruction)\n",
        "      else:\n",
        "        self.recon_loss = self.analytic_recon(self.input_tensor)\n",
        "\n",
        "      # KL Divergence to Gaussian prior in latent space\n",
        "      # We assume prior mean of 0.0 and standard deviation of 1.0\n",
        "      kl = self.gaussian_kl_divergence(self.z_mean, self.expanded_z_logvar,\n",
        "                                       0.0,\n",
        "                                       2 * tf.log(1.0))\n",
        "      avg_kl = tf.reduce_mean(kl, 0)\n",
        "      self.kl_div = tf.reduce_sum(avg_kl)\n",
        "      self.elbo = self.recon_loss + self.kl_div\n",
        "\n",
        "  def analytic_recon(self, x):\n",
        "    \"\"\"Compute the analytic reconstruction error for a linear VAE.\n",
        "\n",
        "    Args:\n",
        "      x: The input tensor.\n",
        "\n",
        "    Returns:\n",
        "\n",
        "      E_q[log p(x|z)]\n",
        "    \"\"\"\n",
        "    wv = tf.matmul(self.enc_weight, self.dec_weight)\n",
        "    x_sub_mu = x - self.mean\n",
        "\n",
        "    wvx = tf.matmul(x_sub_mu, wv)\n",
        "    xvwwvx = tf.reduce_sum(wvx * wvx, 1)\n",
        "\n",
        "    tr_wdw = tf.trace(\n",
        "        tf.matmul(\n",
        "            self.dec_weight,\n",
        "            tf.expand_dims(tf.exp(self.z_logvar), 1) * self.dec_weight,\n",
        "            transpose_a=True))\n",
        "\n",
        "    xwvx = tf.reduce_sum(wvx * x_sub_mu, 1)\n",
        "\n",
        "    xx = tf.reduce_sum(x_sub_mu * x_sub_mu, 1)\n",
        "\n",
        "    d = self._config.input_dim\n",
        "    recon_loss = 0.5 * (\n",
        "        (tr_wdw + xvwwvx - 2.0 * xwvx + xx) / self.sigma_sq + d *\n",
        "        (_LOG_2PI + self.log_sigma_sq))\n",
        "    return tf.reduce_mean(recon_loss)\n",
        "  \n",
        "  def stochastic_recon(self, x, logits):\n",
        "    \"\"\"Stochastic estimation of the reconstruction error.\n",
        "\n",
        "    Computes the reconstruction error under a Gaussian observation model.\n",
        "\n",
        "    Args:\n",
        "      x: The input tensor.\n",
        "      logits: The stochastic output of the decoder.\n",
        "\n",
        "    Returns:\n",
        "\n",
        "      MEAN_i(log p(x_i|z_i))\n",
        "    \"\"\"\n",
        "    input_dim = tf.cast(tf.shape(x)[1], tf.float32)\n",
        "    batch_size = tf.shape(x)[0]\n",
        "    inputs = tf.reshape(x, (batch_size, -1))\n",
        "    logits = tf.reshape(logits, (batch_size, -1))\n",
        "\n",
        "    return tf.reduce_mean(0.5 * (tf.reduce_sum(\n",
        "        (inputs - logits)**2 / self.sigma_sq, 1) +\n",
        "                                  (_LOG_2PI + self.log_sigma_sq) * input_dim))\n",
        "\n",
        "  def gaussian_kl_divergence(self, q_mean, q_logvar, p_mean, p_logvar):\n",
        "    \"\"\"KL Divergence between two Gaussian distributions.\n",
        "\n",
        "    Given q ~ N(mu_1, sigma^2_1) and p ~ N(mu_2, sigma^2_2), this function\n",
        "    returns,\n",
        "\n",
        "    KL(q||p) = log (sigma^2_2 / sigma^2_1) +\n",
        "      (sigma^2_1 + (mu_1 - mu_2)^2) / (2 sigma^2_2) - 0.5\n",
        "\n",
        "    Args:\n",
        "      q_mean: Mean of proposal distribution.\n",
        "      q_logvar: Log-variance of proposal distribution\n",
        "      p_mean: Mean of prior distribution.\n",
        "      p_logvar: Log-variance of prior distribution\n",
        "\n",
        "    Returns:\n",
        "      The KL divergence between q and p ( KL(q||p) ).\n",
        "    \"\"\"\n",
        "    return 0.5 * (\n",
        "        p_logvar - q_logvar + (tf.exp(q_logvar) +\n",
        "                               (q_mean - p_mean)**2) / tf.exp(p_logvar) - 1)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "s6WkR3WIRTKO",
        "colab_type": "text"
      },
      "source": [
        "## Probabilistic PCA"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7zmfyJrWPGl_",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def compute_ppca_mle(covariance, z_dim):\n",
        "  \"\"\"Compute the probabilistic pPCA MLE.\n",
        "\n",
        "  Args:\n",
        "    covariance: The data covariance matrix\n",
        "    z_dim: The number of hidden dimensions\n",
        "  \n",
        "  Returns:\n",
        "    (sigma_sq_mle, W_mle): The MLE solution to pPCA.\n",
        "  \"\"\"\n",
        "  w, u = np.linalg.eigh(covariance)\n",
        "  eigvals, eigvecs = w[::-1], u[:,::-1]\n",
        "  missing_eigvals = eigvals[z_dim:]\n",
        "  sigma_sq_mle = missing_eigvals.sum() / (eigvals.shape[0] - z_dim)\n",
        "  \n",
        "  active_eigvals = np.diag(eigvals[:z_dim])\n",
        "  active_components = eigvecs[:,:z_dim]\n",
        "  \n",
        "  W_mle = active_components.dot(\n",
        "      (active_eigvals - sigma_sq_mle * np.eye(z_dim)) ** 0.5)\n",
        "  return sigma_sq_mle, W_mle\n",
        "\n",
        "def log_p_x_true(W, sigma_sq, N, data_cov):\n",
        "  \"\"\"Compute the model likelihood of pPCA.\n",
        "\n",
        "  Args:\n",
        "    W: The decoder weight.\n",
        "    sigma_sq: Sigma^2 of observation model\n",
        "    N: Total number of data points\n",
        "    data_cov: The dxd data covariance matrix\n",
        "  \"\"\"\n",
        "  d = data_cov.shape[0]\n",
        "  C = W.dot(W.T) + sigma_sq * np.eye(W.shape[0])\n",
        "  loglik = d * _LOG_2PI + np.linalg.slogdet(C)[1]\n",
        "  loglik += np.trace(np.linalg.inv(C).dot(data_cov))\n",
        "  return -loglik * N / 2"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "usQXih68Skcl",
        "colab_type": "text"
      },
      "source": [
        "# Training the Linear VAE and pPCA"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "M56ZOj9YVpXO",
        "colab_type": "text"
      },
      "source": [
        "## Probabilistic PCA"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "a_EelWqOSnc8",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Some hyperparameters\n",
        "DATA_SIZE = 1000\n",
        "Z_DIM = 200\n",
        "VAE_STEPS = 10000"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dYdX2rF7SEUV",
        "colab_type": "code",
        "outputId": "f400288d-16ee-4b7d-dabb-b4dde09f8621",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        }
      },
      "source": [
        "mnist_data = np.array(list(tfds.as_numpy(load_mnist_data(DATA_SIZE))))\n",
        "covariance = np.cov(mnist_data, rowvar=False)\n",
        "\n",
        "sigma_mle, W_mle = compute_ppca_mle(covariance, Z_DIM)\n",
        "avg_log_p_mle = log_p_x_true(W_mle, sigma_mle, DATA_SIZE, covariance)/DATA_SIZE\n",
        "print(\"Avg. Maximum Likelihood: {}\".format(avg_log_p_mle))"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Avg. Maximum Likelihood: 941.6543285914813\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sIWlFb9fVu59",
        "colab_type": "text"
      },
      "source": [
        "## Linear VAE"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Z0yp9prMVEwL",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def build_train_op(mnist_data, learning_rate=0.001):\n",
        "  \"\"\"Builds and returns a training op for optimizing the linear VAE\n",
        "\n",
        "  Resets the default graph and builds a new VAE and optimizer train op.\n",
        "\n",
        "  Args:\n",
        "    learning_rate: Learning rate for the Adam optimizer (default 0.01)\n",
        "  \n",
        "  Returns:\n",
        "    (vae, train_op)\n",
        "  \"\"\"\n",
        "  tf.reset_default_graph()\n",
        "\n",
        "  input_data = tf.data.Dataset.from_tensor_slices(mnist_data)\n",
        "  input_data = input_data.repeat().batch(DATA_SIZE)\n",
        "  input_tensor = input_data.make_one_shot_iterator().get_next()\n",
        "\n",
        "  # Build linear VAE\n",
        "  vae_config = LinearVAEConfig(hidden_dim=Z_DIM, sigma_init=1.0)\n",
        "  vae = LinearVAE(input_tensor, vae_config)\n",
        "\n",
        "  # Make optimizer train op\n",
        "  optimizer = tf.train.AdamOptimizer(learning_rate=learning_rate)\n",
        "  return vae, optimizer.minimize(vae.elbo)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fh3kEEX-VuLb",
        "colab_type": "code",
        "outputId": "e8d6488b-95d2-4bcf-d6c5-dcd81c760741",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 261
        }
      },
      "source": [
        "# Use the same data as pPCA\n",
        "vae, train_op = build_train_op(mnist_data, 0.001)\n",
        "sess = tf.InteractiveSession()\n",
        "sess.run(tf.global_variables_initializer())"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "WARNING:tensorflow:From <ipython-input-7-ecc2ffb1a889>:16: DatasetV1.make_one_shot_iterator (from tensorflow.python.data.ops.dataset_ops) is deprecated and will be removed in a future version.\n",
            "Instructions for updating:\n",
            "Use `for ... in dataset:` to iterate over a dataset. If using `tf.estimator`, return the `Dataset` object directly from your input function. As a last resort, you can use `tf.compat.v1.data.make_one_shot_iterator(dataset)`.\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "stream",
          "text": [
            "WARNING:tensorflow:From <ipython-input-7-ecc2ffb1a889>:16: DatasetV1.make_one_shot_iterator (from tensorflow.python.data.ops.dataset_ops) is deprecated and will be removed in a future version.\n",
            "Instructions for updating:\n",
            "Use `for ... in dataset:` to iterate over a dataset. If using `tf.estimator`, return the `Dataset` object directly from your input function. As a last resort, you can use `tf.compat.v1.data.make_one_shot_iterator(dataset)`.\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "stream",
          "text": [
            "WARNING:tensorflow:From /usr/local/lib/python3.6/dist-packages/tensorflow_core/python/ops/math_grad.py:1375: where (from tensorflow.python.ops.array_ops) is deprecated and will be removed in a future version.\n",
            "Instructions for updating:\n",
            "Use tf.where in 2.0, which has the same broadcast rule as np.where\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "stream",
          "text": [
            "WARNING:tensorflow:From /usr/local/lib/python3.6/dist-packages/tensorflow_core/python/ops/math_grad.py:1375: where (from tensorflow.python.ops.array_ops) is deprecated and will be removed in a future version.\n",
            "Instructions for updating:\n",
            "Use tf.where in 2.0, which has the same broadcast rule as np.where\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jN66TT_4XehM",
        "colab_type": "code",
        "outputId": "a2408476-e635-494f-c82d-70039df3fb55",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        }
      },
      "source": [
        "elbos = []\n",
        "# Note that floating point errors in the analytic ELBO computation\n",
        "# lead to values which may (slightly) exceed the true log likelihood here\n",
        "for i in range(12000):\n",
        "  elbo = -sess.run([train_op, vae.elbo])[-1]\n",
        "  if i % 200 == 0:\n",
        "    elbos.append(elbo)\n",
        "    print(\"ELBO at step {} = {}\".format(i, elbo))"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "ELBO at step 0 = -968.6621704101562\n",
            "ELBO at step 200 = -668.5194702148438\n",
            "ELBO at step 400 = -590.5942993164062\n",
            "ELBO at step 600 = -514.5360107421875\n",
            "ELBO at step 800 = -441.3041076660156\n",
            "ELBO at step 1000 = -369.9105529785156\n",
            "ELBO at step 1200 = -299.8696594238281\n",
            "ELBO at step 1400 = -231.04783630371094\n",
            "ELBO at step 1600 = -163.48062133789062\n",
            "ELBO at step 1800 = -97.22725677490234\n",
            "ELBO at step 2000 = -32.38248062133789\n",
            "ELBO at step 2200 = 31.040729522705078\n",
            "ELBO at step 2400 = 92.99249267578125\n",
            "ELBO at step 2600 = 153.3643341064453\n",
            "ELBO at step 2800 = 212.14169311523438\n",
            "ELBO at step 3000 = 269.2225341796875\n",
            "ELBO at step 3200 = 324.53912353515625\n",
            "ELBO at step 3400 = 377.9089050292969\n",
            "ELBO at step 3600 = 429.6019287109375\n",
            "ELBO at step 3800 = 479.26446533203125\n",
            "ELBO at step 4000 = 526.95068359375\n",
            "ELBO at step 4200 = 572.6626586914062\n",
            "ELBO at step 4400 = 616.456787109375\n",
            "ELBO at step 4600 = 658.244384765625\n",
            "ELBO at step 4800 = 698.238525390625\n",
            "ELBO at step 5000 = 736.5476684570312\n",
            "ELBO at step 5200 = 772.7591552734375\n",
            "ELBO at step 5400 = 806.6507568359375\n",
            "ELBO at step 5600 = 837.297119140625\n",
            "ELBO at step 5800 = 863.6258544921875\n",
            "ELBO at step 6000 = 885.6294555664062\n",
            "ELBO at step 6200 = 903.287353515625\n",
            "ELBO at step 6400 = 916.8137817382812\n",
            "ELBO at step 6600 = 926.583984375\n",
            "ELBO at step 6800 = 933.1580810546875\n",
            "ELBO at step 7000 = 937.2884521484375\n",
            "ELBO at step 7200 = 939.6589965820312\n",
            "ELBO at step 7400 = 940.8868408203125\n",
            "ELBO at step 7600 = 941.466064453125\n",
            "ELBO at step 7800 = 941.73388671875\n",
            "ELBO at step 8000 = 941.84033203125\n",
            "ELBO at step 8200 = 941.8969116210938\n",
            "ELBO at step 8400 = 941.9080200195312\n",
            "ELBO at step 8600 = 941.92138671875\n",
            "ELBO at step 8800 = 941.9198608398438\n",
            "ELBO at step 9000 = 941.9194946289062\n",
            "ELBO at step 9200 = 941.8944091796875\n",
            "ELBO at step 9400 = 941.9273681640625\n",
            "ELBO at step 9600 = 941.9351806640625\n",
            "ELBO at step 9800 = 941.93115234375\n",
            "ELBO at step 10000 = 941.921875\n",
            "ELBO at step 10200 = 941.9248046875\n",
            "ELBO at step 10400 = 941.928466796875\n",
            "ELBO at step 10600 = 941.9375\n",
            "ELBO at step 10800 = 941.9421997070312\n",
            "ELBO at step 11000 = 941.936279296875\n",
            "ELBO at step 11200 = 941.9503173828125\n",
            "ELBO at step 11400 = 941.9246826171875\n",
            "ELBO at step 11600 = 941.950439453125\n",
            "ELBO at step 11800 = 941.936279296875\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "puVqrNdsXxtP",
        "colab_type": "code",
        "outputId": "9c8b36c6-037b-40e6-e0dd-6a16630fcd46",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 298
        }
      },
      "source": [
        "plt.style.use('seaborn')\n",
        "fig, ax = plt.subplots(figsize=(8,4))\n",
        "ax.plot(100*np.arange(len(elbos)), elbos, color='blue', label='Linear VAE')\n",
        "ax.axhline(avg_log_p_mle, ls='--', color='red', label='pPCA MLE')\n",
        "ax.set_xlabel('Training steps')\n",
        "ax.set_ylabel('Training ELBO')\n",
        "plt.legend(frameon=True)"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x7fea16108390>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 10
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAEGCAYAAACErvdRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeVxU5f7A8c+ZjUVRQUHTsnIj9zX3\nfbuW5i8X0EzJe81KRa1cI0u7pZaZ19xyLZfMBaqbFaFpbhliSpmaZpq5CwMhsg3DzJzfH3OdpAFB\nGxhgvu/XixfMmWfO+c5X8DvPeZ7zHEVVVRUhhBBCeASNuwMQQgghRPGRwi+EEEJ4ECn8QgghhAeR\nwi+EEEJ4ECn8QgghhAfRuTuA4mA0prl0f/7+vqSkZLp0n6Wd5CQ3yYczyYkzyYkzyYmzu8lJYKBf\nvs+5pcd/+vRpevbsyYcffgjA1atXGTFiBMOGDWPixImYzWYAtm3bxqBBgwgJCSEyMhKAnJwcJk2a\nxBNPPMHw4cO5ePFiscev02mL/ZglneQkN8mHM8mJM8mJM8mJM1fnpNgLf2ZmJq+//jrt2rVzbFu0\naBHDhg3jo48+4v777ycqKorMzEyWLl3K2rVr2bBhA+vWreP69et88cUXVKhQgU2bNvHcc8/xzjvv\nFPdbEEIIIUqtYi/8BoOBVatWERQU5NgWFxdHjx49AOjWrRuxsbEcPXqUxo0b4+fnh7e3Ny1atCA+\nPp7Y2Fh69eoFQPv27YmPjy/utyCEEEKUWsU+xq/T6dDpch82KysLg8EAQOXKlTEajSQlJREQEOBo\nExAQ4LRdo9GgKApms9nx+rz4+/u6/FTJ7cZPPJXkJDfJhzPJiTPJiTPJiTNX5qTETe7LbwXhO91+\nK1dPFAkM9HP5hMHSTnKSm+TDmeTEmeTEmeTE2d3kpMRN7vsrX19fTCYTAAkJCQQFBREUFERSUpKj\nTWJiomO70WgE7BP9VFW9bW9fCCGEEH8qEYW/ffv2bN++HYAdO3bQqVMnmjZtyrFjx7hx4wYZGRnE\nx8fTqlUrOnToQExMDAC7d++mTZs27gxdCCGEKFWK/VT/8ePHeeutt7h8+TI6nY7t27czf/58pk+f\nzpYtW6hevTqPP/44er2eSZMmMWrUKBRFYdy4cfj5+fHoo4/y3Xff8cQTT2AwGHjzzTeL+y0IIYQQ\npZbiCbfltda8P8/tmWMnYBr1DAB+Y0ejj4t1apPTshVpK9cC4L1hLb4L56PVKFhtudP2R2w8GAxo\nfz1NxaED8zxe2oLF5HTpBkClf3RFc8tQxk2m0CfInPYyAOVmvozXF5/l+X5SP/0SAMNXX1J+xrQ8\nj3f98+3YqtdAuZ6Cf49OebbJiHiV7EGhAFR4MgTdqZNObczdepI+fyEAPosX4rN2tVMbrV95jHsO\nAqA7fIgKz/4rz+PdeH8DlqbNAfBv0wzFYnFqk/XMGLKeHQdA+efHYdi/16mNpXFTbqzdCIDX5o2U\ne3tunsf7Y+9BKF8eze/nqDTosTzbpM9bgLlHbwAq9euN5uoVpzbZAwaTMWMWAOXemIXXp1FObWz3\nVOf6FzsACDz8LdZnn8vzeNc//hzbAw9CejoBXdrm2SZjyktkD30SgAojn0R37KhTG3OnLqQvXAqA\nz4ql+Kx8z6mNqtOREvcjALqjP1DhXyPyPN6NFe9jadUaAP9OrVEynefFZI18mqzxzwNQfvLzGHbv\ndGpjeag+Nzba19zw+ngr5eb82/HcrX83Kbv2o1byR3PlMpUe+0eeMaW/8RbmR/oCUHFAX7QXzju1\nye73f2S8NhsA37dm4711k1MbW5UqXN++B7MZLDF7qDYjHJsK6l++to35jAT/h7Bkmnnmncagggq5\nvm9vNY1v6z+NzQb/+no4DyYcsh9E/fPb6aAOrOj4AaoKPU69R9+T/8kd0P/ajn/sV7y99QRei2fq\n/oF5NWFZm/c5GWj/231rRyt8zalO7293rX/ySYMIAEbGv0DLK184tbnqV485Xez/Z7S5+DHDj053\nagMwo+e3pHpXpaIpgTd2dsyzzYdN3yTuvkEAROztyz1pp53aHKnej7Ut7O974M9z6PbbB05tMg0V\nmdb7MAD1jfsZG2f/P0NR7P8eN83r+AkXKzUGYPEXdfOM6fPgF9lRdwwA4+JG8pDxgFObM5Vb8247\n+/8Z3c+uYcDJvDuNzz9yAqvWwD2pp4jYn/f/GStaLedYVfuVaHN2tqdCttGpzd4HhhPZcCYAI45O\no82lT5zaJJZ7gH93/Rq9XuXlRp8S/N7UPI+n/e4ARq+Kd/R/udeOmDzbQQmc3CeEEIWhqpCeDj/E\na7h2TUOdH7XUTwWrVcFqBZvN/mW8qqXD/eXJylLogQ+rUfLc3+tvePML3ujREJpPm29261m92z6n\nqAsaAvNo89s5DRvO2dtUQEfXfOLfskUPQCP0jMunzTe7dezH3u6lfGI6dlzL1uP2Nh3QUD+PNtcy\nFbZutbexoiPvrgl8tk1PInqC0BORT5vvYnV8HGvf179QqJhHm1/PaNh6xt6mLlryGpC9fktMndAx\nMp/jbd+h5/j/cvBWPm3if9Cy9Qd7m35oqJlHmwuZGrZetLepgJa8P25C1Md6ctATjJ7n82mzd5+O\nXf+LaTJKnoX055+1RP5sb9MaDY3zaGPMVIiMtLcZFKYhOJ/juZpH9PhdPUNUZp06k5zkJvlwdjc5\nsdng8mWFs2c1/Pab/evmzxcuKFiteRdDAK1WpVIlFT8/qFhRpUIF1fHd1xe8vVW8vcHbG3x97T97\neakYDKDVgl6votff/Nn+XatV0WjsP2s09t6pVmv/nvtLdfwM+X+vUqU8ycnpecav5P/WipWqFm8s\nlSvnn5Pidqfv+25zpdNBtWr5l2JXz+qXHr8QokRQVTh3TuHHH7X8+KOWo0c1/PSTlowM5/9Jq1Sx\n0aKFjfvusxEUpFK1qo1q1VSqVr35ZcPPr+QUz/wEBkK5cmW+73VHAgPtH8RE0ZHCL4RwC7MZjhzR\nsnevlu+/1/LTT1pSU/+s1IqiUq+ejQYNbNSqZaN2bfv3WrVsVKrkxsCFKOWk8AshioXNBidOaNi3\nT8u+fTpiY7VkZv5Z6GvXttGjh4Vmzaw0a2ajUSMr5cu7MWAhyigp/EKIImM2w759Wj7/XM+uXZCY\nWM7xXL16Vjp3ttK5s4W2ba3SixeimEjhF0K4VHY27N2rZds2PTExOm7csPfqq1WDkJAcOne20Lmz\nlXvukXFcIdxBCr8Q4m+z2ew9+8hIe7FPS7MX++rVbTzxRA6PPZbDI4+UIznZ5OZIhRBS+IUQdy0h\nQWHTJj0ffqjnwgX7CuD33mvjySdz6N8/hxYtbGj+tzC4pkQsEC6EkMIvhLgjVivs3q1lwwY9O3bo\nsFoVfH1VnngihyefNPPww7YSfxmdEJ5MCr8QolDS0mD9ej1r1hi4dMnefW/c2MqIETkMHJhDhQpu\nDlAIUShS+IUQt3XtmsKqVXrWrjWQlmbv3Y8YYSYsLIemTW3uDk8IcYek8Ash8vTrrxqWLdMTGanH\nbFYIDLQxYYKZp54yy6V3QpRiUviFELkcParhnXcMxMTYbx5Su7aNsWOzCQnJwdvbzcEJIf42KfxC\nCACOH9fw9tsGvvrKXvBbtrQSHm6mTx8LWq2bgxNCuIwUfiE83KlT9oL/+ef2gt+qlZXp07Pp1Mkq\ns/OFKIOk8Avhoc6cUZg/34tPP9WhqgrNmtkLfrduUvCFKMtKROGPjIxk27ZtjsfHjx+nUaNGZGZm\n4uvrC8C0adNo1KgRq1evJiYmBkVRCA8Pp0uXLu4KW4hSKSFBYd48Axs36rHZFBo1sjJtWja9e0vB\nF8ITlIjCHxISQkhICACHDh3iq6++4syZM8ydO5d69eo52l28eJHo6Gg2b95Meno6w4YNo2PHjmhl\nAFKIAqWnw5IlBpYvN5CZqVC3rpXp08307WuRVfWE8CAl7s996dKljB07Ns/n4uLi6NSpEwaDgYCA\nAGrUqMGZM2eKOUIhShezGdas0dO6dTkWLPCifHmV+fNN7N2byWOPSdEXwtOUiB7/TT/99BP33HMP\ngYGBACxatIiUlBRq165NREQESUlJBAQEONoHBARgNBoJDg6+7X79/X3R6Vx7ViAw0M+l+ysLJCe5\nuTsfqgoffwwvvQRnzkD58vDvf8OLL2ooV84bKP5r89ydk5JIcuJMcuLMlTkpUYU/KiqKAQMGABAW\nFkZwcDA1a9Zk5syZbNy40am9qhbutp4pKZkujTMw0A+jMc2l+yztJCe5uTsfx45pePllLw4e1KHT\nqYwalcOLL5oJDFTJzIRM1/5JFIq7c1ISSU6cSU6c3U1ObvdBoUSd5IuLi6N58+YA9OrVi5o1awLQ\nvXt3Tp8+TVBQEElJSY72CQkJBAUFuSVWIUqixESFF1/0omdPXw4e1NGnTw7ffpvB3LnZBAYW7oOy\nEKJsKzGFPyEhgXLlymEwGFBVlZEjR3Ljxg3A/oGgbt26tG3blj179mA2m0lISCAxMZE6deq4OXIh\n3M9shqVL9bRtW44PPzQQHGwjMjKT9etN1KolBV8I8acSc6rfaDQ6xu8VRSE0NJSRI0fi4+ND1apV\nGT9+PD4+PoSGhjJ8+HAURWHWrFloZGaS8HA7dmh59VVvfvtNg7+/yptvmggLy0FXYv66hRAliaIW\ndqC8FHP1eJGMQTmTnORWHPn4/XeFGTO82bFDh1ar8q9/5TB5cjb+/kV62LsmvyPOJCfOJCfOXD3G\nL30CIUoZk8l+Pf6iRQZMJoWOHS3MmZPNQw/JLXKFEAWTwi9EKbJzp5aXXvLm/HkNVavaePddE48/\nbpEV94QQhSaFX4hS4MIFhRkzvIiJ0aPVqowZY2bKlGzKl3d3ZEKI0kYKvxAlmMUCy5frefttL7Ky\nFNq3tzB3bjb168tpfSHE3ZHCL0QJFR+vYdIkb06c0FKlio35800MHiyn9YUQf48UfiFKmLQ0mDvX\nizVr9KiqwrBhZmbOLLmz9YUQpYsUfiFKkOhoHS+95MXVqxrq1LEyf3427dtb3R2WEKIMkcIvRAmQ\nkKAwfboXX36px2BQmTw5m4kTzXh5uTsyIURZI4VfCDdSVYiM1DFjhjfXryu0bWvhnXeyqVtXJu8J\nIYqGFH4h3OTKFYUpU7z5+msdvr4qc+ea+Oc/c5BVqIUQRUkKvxDFTFXho4/0vPqqF2lpCp07W1iw\nwETNmmV+9WwhRAkghV+IYnTpksILL3izd6+O8uVV3nnHxPDhOXKJnhCi2EjhF6IYqCps2aIjIsKb\n9HSFHj0szJ9vokYN6eULIYqXFH4hipjRqDB5shdffaWnfHmVhQuzeOIJWYhHCOEeUviFKEJffaVj\n0iQvkpI0tG9vYdEiGcsXQriXFH4hikBqKkyY4M3mzXq8vFRef93E6NEyY18I4X5S+IVwsdhYLePH\nw4ULepo0sbJ0qYngYLkuXwhRMpSIwh8XF8fEiROpW7cuAPXq1ePpp59m6tSpWK1WAgMDefvttzEY\nDGzbto1169ah0WgIDQ0lJCTEzdELYZeTA/PnG1i40IBGA5MmZfPii2b0endHJoQQfyoRhR+gdevW\nLFq0yPH4pZdeYtiwYTzyyCMsWLCAqKgoHn/8cZYuXUpUVBR6vZ7BgwfTq1cvKlWq5MbIhYBz5xTG\njvXhyBEtNWva2LxZoU4ds7vDEkIIJyV2xDEuLo4ePXoA0K1bN2JjYzl69CiNGzfGz88Pb29vWrRo\nQXx8vJsjFZ5MVWHrVh3du5fjyBEtgwfnsHt3Bu3auTsyIYTIW4np8Z85c4bnnnuO1NRUwsPDycrK\nwmAwAFC5cmWMRiNJSUkEBAQ4XhMQEIDRaCxw3/7+vuh0WpfGGxjo59L9lQWelpPUVBgzBjZtAj8/\n2LABhg/XA/Zz+56Wj8KQnDiTnDiTnDhzZU5KROF/4IEHCA8P55FHHuHixYuEhYVhtf55K1JVzfvy\np/y2/1VKSqZL4rwpMNAPozHNpfss7TwtJ4cPa3juOR8uXNDQsqWV997L4oEHVG5+DvW0fBSG5MSZ\n5MSZ5MTZ3eTkdh8USsSp/qpVq/Loo4+iKAo1a9akSpUqpKamYjKZAEhISCAoKIigoCCSkpIcr0tM\nTCQoKMhdYQsPZLPB0qV6+vf35eJFhRdfzGbbtkweeECuzRdClA4lovBv27aNNWvWAGA0GklOTmbg\nwIFs374dgB07dtCpUyeaNm3KsWPHuHHjBhkZGcTHx9OqVSt3hi48yB9/QFiYD6+95k1AgEpUVBbT\np8usfSFE6VIiTvV3796dyZMns2vXLnJycpg1axb169dn2rRpbNmyherVq/P444+j1+uZNGkSo0aN\nQlEUxo0bh5+fjAWJonfokIZnn/Xh8mUNnTtbWLbMRFCQ9PKFEKWPohZ2oLwUc/V4kYxBOSurObGf\n2jcwZ44BVYUpU8w8/7wZbQFzRctqPv4OyYkzyYkzyYkzV4/xl4gevxAl0R9/QHi4Dzt36qha1caK\nFSbat7cW/EIhhCjBpPALkYcff9QwapQPFy9q6NrVwtKlJgIDy/zJMSGEBygRk/uEKEk+/FBPv36+\nXLqkMGVKNps3Z0nRF0KUGdLjF+J/srLgpZe8+OgjA5UqqSxfnkX37nJqXwhRtkjhFwL4/XeFUaN8\nOHZMS9OmVtasyaJmTenlCyHKHjnVLzzezp1aevcux7FjWoYPN/P555lS9IUQZZb0+IXHstlg4UID\nb71lwGCAhQuzGDbM4u6whBCiSEnhFx4pPR3Cw72JjtZz77021q7NokkTm7vDEkKIIieFX3ic335T\neOopH375RUuHDhZWrTJRpYqc2hdCeAYZ4xceZdcu+3j+L79oeeYZM1u3ZknRF0J4FOnxC4+gqrB4\nsYHZs+3j+YsXZzFkiIznCyE8jxR+UeZlZsLEid589pme6tXt4/nNmsl4vhDCM0nhF2XalSsKYWE+\n/PSTlrZtLaxeLXfVE0J4ttsW/t9++4333nuPkydPoigKTZo0YcyYMdx7773FFZ8Qdy0+XkNYmA+J\niRqGDzfz5pvZGAzujkoIIdwr38l9sbGxjBo1iiZNmjB37lzmzJlDw4YNGTVqFPHx8cUZoxB37JNP\ndDz+uC9JSQqvv27inXek6AshBNymx79s2TJWrFhBvXr1HNsaN25MmzZt+Pe//826deuKJUAh7oTN\nBvPmGViwwAs/P5UPPsiiRw9Zb18IIW7Kt/CbTKZcRf+m2rVrk5WVVaRBCXE3MjLsi/J8+aWe+++3\nsXFjFvXqySQ+IYS4Vb6FPycnJ98X3e65uzVv3jyOHDmCxWLh2Wef5ZtvvuHEiRNUqlQJgFGjRtG1\na1e2bdvGunXr0Gg0hIaGEhIS4vJYROlz7ZrCk0/ab7LToYOFNWuyCAhwd1RCCFHy5Fv477vvPnbu\n3EnPnj1zbY+OjubBBx90aRAHDx7k119/ZcuWLaSkpDBgwADatm3Liy++SLdu3RztMjMzWbp0KVFR\nUej1egYPHkyvXr0cHw6EZzp+XMPw4T5cuSKT+IQQoiD5Fv4pU6bw9NNPEx0dTZMmTbDZbPzwww+c\nOXOGjRs3ujSIhx9+mCZNmgBQoUIFsrKysFqdx2WPHj1K48aN8fPzA6BFixbEx8fTvXt3l8YjSo+d\nO7WMHu1DRobCq6+aGDcuB0Vxd1RCCFFyKaqq5ntRs9ls5r///S8///wzPj4+BAcH07dvX/R6fZEF\ntGXLFg4fPoxWq8VoNJKTk0PlypV55ZVXOHDgAMeOHSMiIgKAhQsXcs899zBkyJDb7tNisaLTaYss\nZuEey5bB+PFgMMCGDTB4sLsjEkKIku+21/EbDAZCQ0MByM7O5uzZs2RnZxdZ4d+5cydRUVG8//77\nHD9+nEqVKlG/fn1WrlzJkiVLaN68ea72t/nMkktKSqZL4wwM9MNoTHPpPku74syJ1QqzZnmxYoWB\nKlVsbNiQRcuWNozGYjl8ocjviDPJiTPJiTPJibO7yUlgoF++z+V7Hf+pU6cYN24cc+bM4erVq/zf\n//0fL7/8Mr1792bPnj13FEBh7N+/n+XLl7Nq1Sr8/Pxo164d9evXB6B79+6cPn2aoKAgkpKSHK9J\nTEwkKCjI5bGIkisjA/75T29WrDBQr56Vr77KpGVLmbkvhBCFlW/hf+ONN+jZsycVK1Zk1KhRvPPO\nO3z66ad8+umnLF++3KVBpKWlMW/ePFasWOGYqDd+/HguXrwIQFxcHHXr1qVp06YcO3aMGzdukJGR\nQXx8PK1atXJpLKLkSkxUGDDAl5gYPZ06Wfjyy0zuv1+W3xVCiDuR76l+RVEYMGAAADExMTRs2BCA\nqlWruvxUf3R0NCkpKTz//POObQMHDuT555/Hx8cHX19f5s6di7e3N5MmTWLUqFEoisK4ceMcE/1E\n2Xb2rMKQIb5cuKBh6NAc5s83ycx9IYS4C7ct/Df5+/vn+5wrDBkyJM8Jejc/eNyqT58+9OnTx6XH\nFyXb999rGDHChz/+0DB5cjZTpphl5r4QQtylfAt/YmIiUVFRABiNRsfPNx8LURyio3U895w3OTmw\nYIGJ4cNdv3iUEEJ4knwLf/PmzTly5AgAzZo1c/x887EQRW3NGj0vv+yFtzds2JBFz56y5r4QQvxd\n+Rb+uXPn5vuiGzduFEkwQoD9RjuzZxtYvNiLKlVsfPRRFs2aycx9IYRwhXxn9d9OeHi4q+MQAoCc\nHPuNdhYv9qJWLRvR0ZlS9IUQwoVuu4BPfgq7cI4QdyI9HUaN8mH3bh0tW1r58MMsKleW3zUhhHCl\nuyr8rp7VL0Rysv3uevHxWnr2tLBqVRblyrk7KiGEKHvyLfyxsbH5vkjG+IUrXbyoMGSID2fOaAkN\nzeE//zFRhLeDEEIIj5Zv4V+2bFm+L5JFc4SrnDypYehQH65e1TBunJlXX82Wa/SFEKII5Vv4N2zY\nUJxxCA908KCWESN8SE1VmDXLxNixco2+EEIUtXxn9b/11lu5Hu/evdvx8zPPPFN0EQmPsGOHltBQ\nHzIyYOnSLCn6QghRTPIt/CdOnMj1+IMPPnD8nJWVVXQRiTIvKkrHU0/5oCj2hXlCQizuDkkIITxG\noa/jv/USPpnVL+7WmjV6xo71oVw5iIzMpEcPWY1PCCGKk1zHL4qFqsJ//mPgzTe9CAy0sWVLFo0a\nycI8QghR3PIt/KqqOr7+uk2IO6GqMHOmF8uXG7jvPhuRkZnUqiW/R0II4Q75Fv7vv/+eBg0aOB6r\nqkqDBg1QVVVO9YtCs1hg0iRvNm3SU6+ela1bs6heXYq+EEK4S76F/9SpU8UZhyiDsrPhuee8+fJL\nPc2aWdm0SZbgFUIId7urMX4hCpKZCSNH+rBnj44OHSysX5+FrPskhBDuVyoL/5w5czh69CiKohAR\nEUGTJk3cHZK4RVoaPPmkDwcP6ujVy8Lq1Vn4+Lg7KiGEEHCXt+V1p0OHDnH+/Hm2bNnC7NmzmT17\ntrtDErdISYGQEF8OHtTRv38OH3wgRV8IIUqSAnv8UVFRzi/S6XjwwQdp2rRpkQR1O7GxsfTs2ROA\n2rVrk5qaSnp6OuXLly/2WERuiYkKoaE+/PyzliFD7Dfb0ZXKc0pCCFF2Ffjf8oEDBzhw4AAtWrRA\nq9Vy5MgRHn74YS5evEiXLl144YUXiiNOh6SkJBo2bOh4HBAQgNFovG3h9/f3RafTujSOwEAZsL7V\npUswaFB5fvkFxoyBJUv0aDSefYs9+R1xJjlxJjlxJjlx5sqcFFj4rVYr0dHRVKlSBYDk5GTmzp3L\np59+ytChQ10WyN0qzLoCKSmZLj1mYKAfRmOaS/dZmv3+u0JoaHl+/x3HHfaSk90dlXvJ74gzyYkz\nyYkzyYmzu8nJ7T4oFFj4ExISHEUfoHLlyly6dAlFUbDZin/ltaCgIJKSkhyPExMTCQwMLPY4hN2Z\nMwqDBvly9SpMnZrNpElmua2uEEKUYAUW/urVqzNhwgRat26Noij88MMPlCtXjpiYGO65557iiDGX\nDh06sHjxYoYOHcqJEycICgqS8X03OXVKw6BBPhiNGubPh7Aws7tDEkIIUYACC/9bb73FZ599xqlT\np7DZbDRt2pQBAwaQkZFBly5diiPGXFq0aEHDhg0ZOnQoiqIwc+bMYo9BwPHjGkJCfEhO1jB3rolJ\nk7wxGt0dlRBCiIIUWPgNBgN9+vShbdu2jm0pKSncd999RRrY7UyePNltxxbw008aQkJ8uX4d5s83\nERaWA3i7OywhhBCFUGDhf+ONN/j4448JCAgAcKzVv2vXriIPTpQ8R45oGDLEl7Q0ePddE0OHWtwd\nkhBCiDtQYOGPi4vj4MGDeHl5FUc8ogSLi9PyxBM+ZGXBsmUmBg2Soi+EEKVNgYX//vvvl6IvOHBA\ny5NP+mA2w4oVJvr3l6IvhBClUYGFv1q1ajz55JO0bNkSrfbPRXAmTpxYpIGJkmP/fi3Dh/tgscCa\nNSYeeUSKvhBClFYFFv5KlSrRrl274ohFlED79mkZMcIHqxXWrs2iVy+ru0MSQgjxN+Rb+G9O4hs7\ndmxxxiNKkH377D19m81e9Hv2lKIvhBClXb6F/6mnnmL9+vU0aNAA5Zal2G5+IDh58mSxBCjcY+9e\ne0/fZoN167Lo0UOKvhBClAX5Fv7169cDcOrUqWILRpQMe/ZoCQvzQVWl6AshRFlT4Bi/0WgkOjqa\n1NTUXDfEkcl9ZdPu3VqeeurPot+9uxR9IYQoSzQFNXj22Wc5deoUGo0GrVbr+BJljxR9IYQo+wrs\n8fv6+jJ37tziiEW40d69UvSFEMITFNjjb9q0KWfPni2OWISbfPutfUz/5kQ+KfpCCFF2Fdjj379/\nP2vXrsXf3x+dTueY1b9nz55iCE8UtdjYPxfnkaIvhBBlX4GF/7333iuOOIQb3Fx7PycH3n9frtMX\nQghPkG/h37t3L126dCE2NjbP5wcPHlxkQYmi9/33GoYOta+9v2qViX/8Q4q+EEJ4gnwL/y+//EKX\nLl04cuRIns9L4S+94uM1DB3qi8kEK1ea6NtX1t4XQghPkW/hf+aZZwDynNF/c3EfUfocPaohNNSX\njAxYvtzEY49J0RdCCE9S4MZkB44AACAASURBVBj/yZMnWb58OSkpKQCYzWauXbtGWFiYSwKwWCy8\n/PLLXLhwAavVytSpU2nVqhUjRowgMzMTX19fAKZNm0ajRo1YvXo1MTExKIpCeHg4Xbp0cUkcnuD4\ncQ0hIb6kpcHSpSYef1yKvhBCeJoCC/9rr73GiBEjWLlyJS+88AIxMTG8+OKLLgvgs88+w8fHh02b\nNvHrr7/y0ksvERUVBdjPNtSrV8/R9uLFi0RHR7N582bS09MZNmwYHTt2lAWFCuGXXzSEhvpw/brC\nokVZDB4sRV8IITxRgdfxe3t707dvX/z8/OjatSuzZ89mzZo1Lgugf//+vPTSSwAEBARw/fr1fNvG\nxcXRqVMnDAYDAQEB1KhRgzNnzrgslrLqt98UBg/2ISlJw9tvmxg6VIq+EEJ4qgJ7/NnZ2Zw+fRov\nLy8OHTpEnTp1uHz5sssC0Ov1jp/XrVtHv379HI8XLVpESkoKtWvXJiIigqSkJAICAhzPBwQEYDQa\nCQ4Ovu0x/P190elce1YgMNDPpfsrKr//DiEhkJAACxfCxInegHeRHKu05KS4SD6cSU6cSU6cSU6c\nuTInBRb+yZMnc/HiRSZMmMDUqVNJTk5m9OjRd3WwyMhIIiMjc20bP348nTp1YuPGjZw4cYLly5cD\nEBYWRnBwMDVr1mTmzJls3LjRaX+33jTodlJSMu8q3vwEBvphNKa5dJ9F4coVhf79fbl4UcOMGdkM\nG2bGaCyaY5WWnBQXyYczyYkzyYkzyYmzu8nJ7T4oFFj4fXx8aNmyJQDbt2+/owP/VUhICCEhIU7b\nIyMj+eabb1i2bJnjDECvXr0cz3fv3p3o6GjatGnDuXPnHNsTEhIICgr6WzGVVQkJCoMG+XLhgobJ\nk7OZMMHs7pCEEEKUAAWO8b/55ptFGsDFixfZvHkzS5YswcvLC7D35EeOHMmNGzcA+9h+3bp1adu2\nLXv27MFsNpOQkEBiYiJ16tQp0vhKo+RkhZAQH86e1TB+fDZTpkjRF0IIYVdgj7969eqMGDGCpk2b\n5hqPnzhxoksCiIyM5Pr16451AwDWrFlDaGgoI0eOxMfHh6pVqzJ+/Hh8fHwIDQ1l+PDhKIrCrFmz\n0GgK/OziUVJTITTUh1OntIwebWbGDDOK4u6ohBBClBSKWsBA+ZIlS/LcHh4eXiQBFQVXjxeV1DGo\njAwIDfXl+++1jBhhZv787GIr+iU1J+4i+XAmOXEmOXEmOXFWbGP827Zto3///qWqwHsykwnCwnz4\n/nstAwfmMG9e8RV9IYQQpUe+58lvLqIjSr6cHBg92of9+3X06ZPD4sUmZE0jIYQQeZEB8lLOaoXw\ncG+2b9fRubOFlStN3DIVQwghhMgl31P9P/zwA127dnXarqoqiqKwZ8+eIgxLFIaqwpQpXnz6qZ7W\nrS2sW5eFd9GszSOEEKKMyLfwN2jQgAULFhRnLOIOqCq8+qoXH35ooEkTKx99lEW5cu6OSgghREmX\nb+E3GAzUqFGjOGMRd+Dttw2sWGEgONjKli1ZVKjg7oiEEEKUBvmO8Tdp0qQ44xB3YMUKPfPne3H/\n/TYiI7OoXLlwSxcLIYQQ+Rb+KVOmFGccopA2b9bxyiveVK1qIzIyk2rVpOgLIYQoPJnVX4pER+t4\n/nlvKlVS2bo1iwcekKIvhBDizkjhLyX27dPyzDPeeHvDpk2Z1K9vc3dIQgghSiEp/KXAkSMawsJ8\nAFi3LouWLaXoCyGEuDsF3qRHuNfJkxqGDfPFZILVq0106WJ1d0hCCCFKMSn8Jdj58wqhoT6kpCi8\n+24W/fpZ3B2SEEKIUk5O9ZdQRqNCaKgvCQkaXnvNxBNPSNEXQgjx90nhL4HS02HYMB/OndMwcWI2\nY8bkuDskIYQQZYQU/hImOxtGjvTh6FEtw4aZiYgwuzskIYQQZYjbx/g/+eQT3n33XWrWrAlA+/bt\nGTNmDKdOnWLWrFkABAcH89prrwGwevVqYmJiUBSF8PBwunTp4q7QXc5mg/Hjvdm3z3573fnzs1EU\nd0clhBCiLHF74Qd49NFHmTZtWq5ts2fPJiIigiZNmjBp0iT27t1LrVq1iI6OZvPmzaSnpzNs2DA6\nduyItgzcfF5VYcYML/77Xz1t2lhYscKErkT86wghhChLSuSpfrPZzOXLlx33C+jWrRuxsbHExcXR\nqVMnDAYDAQEB1KhRgzNnzrg5Wtd4910Dq1cbqF/fyoYNWfj4uDsiIYQQZVGJ6FMeOnSIUaNGYbFY\nmDZtGpUrV6bCLbebq1y5MkajkUqVKhEQEODYHhAQgNFoJDg4+Lb79/f3Radz7VmBwEA/l+1r9WqY\nMwdq1oSvv9ZSo4br9l2cXJmTskDy4Uxy4kxy4kxy4syVOSnWwh8ZGUlkZGSubX379mX8+PF07dqV\nH374gWnTprF69epcbVQ17zXp89v+VykpmXcXcD4CA/0wGtNcsq8dO7Q8+6wPAQEqmzZlYTDYMBpd\nsuti5cqclAWSD2eSE2eSE2eSE2d3k5PbfVAo1sIfEhJCSEhIvs83b96cP/74A39/f65fv+7YnpCQ\nQFBQEEFBQZw7d85pe2n1448annnGBy8v2Lgxi7p1ZSleIYQQRcvtY/yrVq3iiy++AOD06dMEBARg\nMBioVasWhw8fBmDHjh106tSJtm3bsmfPHsxmMwkJCSQmJlKnTh13hn/Xzp9XGDbMB5MJVqyQ9feF\nEEIUD7eP8T/22GNMmTKFzZs3Y7FYmD17NgARERG8+uqr2Gw2mjZtSvv27QEIDQ1l+PDhKIrCrFmz\n0Gjc/tnljqWk2BfoSUrSMHeuiT59ZP19IYQQxUNRCztQXoq5erzo74xBZWdDaKgPsbE6xo41M2tW\ntktjcxcZl8tN8uFMcuJMcuJMcuLM1WP8pa+7XIrZbDBhgjexsTr698/h1VfLRtEXQghRekjhL0Zz\n5hj49FM9rVtbWLLERCkcpRBCCFHKSekpJmvX6lm0yItatWysX5+Ft7e7IxJCCOGJpPAXg127tEyf\n7kWVKjY2bcrkljWIhBBCiGIlhb+InTih4emnfTAYYP36LB58sMzPpRRCCFGCuf1yvrIsIUFh+HAf\nMjIU1qzJolUruVZfCCGEe0mPv4hkZMCIET5cvqxhxoxsHnvM4u6QhBBCCCn8RcFmg3HjvPnxRy3D\nhpkZP97s7pCEEEIIQAp/kXj9dS+io/V07Ghh3rxsFMXdEQkhhBB2UvhdbMMGPUuXGqhTx8r772dh\nMLg7IiGEEOJPUvhdaO9eLVOnehEQYGPjxiwqVXJ3REIIIURuUvhdJDUVnn7aB60W1q41yWV7Qggh\nSiQp/C4SF6clNVVhzBgzbdvK3faEEEKUTFL4XeTIES0A7dpJ0RdCCFFyyQI+LnL4sL3wt2wphV8I\nIfJis1mx2W6/kJnZbMZiySmmiEqH/HKi0WjQaLR3vD/p8buA1Qrx8VqCg61UrOjuaIQQouQxmTIL\nVdBTUzOLIZrSJb+cWCw5mEx3ni/p8bvAqVMaMjIU6e0LIUQebDYrGo0Gg6Hg25LqdBpAlje/Vf45\n0WM2m/6X38L3/N1e+N977z2+++47AGw2G0lJSWzfvp3u3btTrVo1tFr7m5k/fz5Vq1Zlzpw5HD16\nFEVRiIiIoEmTJu4MH/jzNL+sxS+EEM5sNttdnZIWBdNotHecX7cX/jFjxjBmzBgAPv30U5KTkx3P\nrVq1inLlyjkeHzp0iPPnz7NlyxbOnj1LREQEW7ZsKfaY/+rPwi89fiGEECVbiRnjt1gsbNq0ieHD\nh+fbJjY2lp49ewJQu3ZtUlNTSU9PL64Q83X4sBY/P5V69aTHL4QQomRze4//ph07dtCxY0e8vf8c\nA5o5cyaXL1+mZcuWTJo0iaSkJBo2bOh4PiAgAKPRSPny5W+7b39/X3Q6155mCgz0AyA5Gc6ehV69\noGpVP5ceo7S5mRNhJ/lwJjlx5gk5MZvNpKZm/m+sumCFbXcnrly5QkTEFNau3Zhr+3/+8zZDhgyj\nevUaLj8mwMGD3xEZuYV33nnXsS0jI52hQwfx6adfoNPpOXbsKKNH/5P16zdRr14wAF98sY2VK9+j\nRo17Ha+rVq0aM2e+/pcjaKhY0RfDHawPX6yFPzIyksjIyFzbxo8fT6dOnfj444957bXXHNsnTJhA\np06dqFixIuPGjWP79u1O+1PVwq2Ol5Li2lmigYF+GI1pAOzcqQV8adIkG6PRc+/Cd2tOhOQjL5IT\nZ56Skz9n8xd8VlSn02CxuP7sqdVqQ1Vx2vf48ZMA5+2u0rz5w8yZ8zopKan4+dk/5O3evZt27ToC\nWiwWGzExX1Gz5v1s3x5DrVp1AbDZVLp370V4+PO5cvLXOC0WG8nJ6eh0+lzbb/eBslgLf0hICCEh\nIU7bMzMzuXbtGvfe++cnm8cff9zxc+fOnTl9+jRBQUEkJSU5ticmJhIYGFi0QRfg5vj+ww/L+L4Q\nQhTGrFlefP65a8vPY49ZmDUr+45fFx7+DC++OJXdu3eRkZHOhQvnuXz5EhMmTKJduw7s3fsNmzd/\niFarIzi4PuPHv0BGRjqvvTaDrKwsTCYTL7wwhQYNGjF06ADatu2Av78/Tz01CgCtVkunTl3Yv38P\njz76GADffLOTJ56wD2tbrVb27PmG116bw+zZsxgzZrzrkpKPEjHGf+rUKWrVquV4nJaWxqhRozCb\n7T3o77//nrp169KhQwdHz//EiRMEBQUVeJq/qN0s/C1aSOEXQojSLDExgfnzFzFx4mS2bfuEzMxM\n1q1bw7vvLmfJkpUkJibw008/kpycTL9+j7N48Qqeey6cjRvXAfa5am3btncU/Zt69erDN998DUB6\nejrnz5+jWbMWABw+fIgHHniQZs1aUKFCRY4f/6nI32eJGOM3Go0EBAQ4Hvv5+dG5c2eGDBmCl5cX\nDRo0oE+fPiiKQsOGDRk6dCiKojBz5kw3Rv3nwj1161rlTnxCCFFIs2Zl59s7L6pT/YXRpEkzAIKC\ngkhPT+fcud9ISLjGiy+GA/ax+WvXrlGrVh3WrVvNpk0byMnJyTU3rUGDhk77bdSoCZcvX+LGjVQO\nHNhP587dUBQFgK+/jqFnz38A0KvXP9i5czuNGtkvU//mm685depnFEVBVVV69OjNgAGD//b7LBGF\n/x//+Af/+Mc/cm176qmneOqpp5zaTp48ubjCKtAvv2hIT1fk+n0hhCgDbq4bA/Y5ZHq9/fT+ggVL\ncrV7//2VVKkSxCuvvM6pUz+zZMlCx3N/HWu/qWvXHuzbt4d9+/bw9NPPApCdnc233+7jl19O8vHH\nW7FYckhLS2PCBPu8g7zG+F2hRJzqL61kfX4hhCi7atZ8gN9/P0dKyh8ArFmzAqMxkdTU647Z9nv3\n7sZisRS4r169+rB79y6SkhKpV+8hAA4c2E/Llq3YsGEra9d+xIcfRnL//Q8QH3+46N4UJaTHX1rd\nvCOfLNwjhBAl34UL5wkPf8bxeOzYCbdt7+3tzcSJk5g8eSIGg566dYOpUiWQPn368sYbM9m9eyeD\nBoWyc+cOvvxy2233VatWbZKTk+jatbtj29dfx9Cv3//lavfoo4+xa9cOGjdu6nSqH+A//1mKXp/3\nWYXCUtTCXhNXirn6cpmbl+B06ODL1asafv01Ha2Hr0bpKZclFZbkw5nkxJmn5OTm5Xz5nQa/lTvH\n+Euq2+Ukv9ze7nI+OdV/l1JS4NdftbRoYfX4oi+EEKL0kMJ/l374QU7zCyGEKH2k8N+l77+Xwi+E\nEKL0kcJ/l2ThHiGEEKWRFP67YLPZF+6pXdvGLesOCSGEECWeFP67cPIkpKUpcppfCCFEqSPX8d+F\n2Fj7dyn8QghRus2ePYtffjlJhQoVUVWVgIDKvPTSK/j6liMrK4tFixbwyy8/YzB4UaFCBSZNmk7V\nqtUcrx82bBBt2rRn4sRJTvuOjz/M88+P5ZNPvqRKFfsN5axWKwMGPMr//d9ARo161nGToFq16jhe\nd/XqFcLChhIcbF/o5+Z1/HPmvE2FChX/9nuWwn8XpPALIUTZ8eyz4XTo0AmwL8e7desmRo58mkWL\nFnDPPfcwbdrLgP2uerNmRfDee+8DcOrUSVRVZc+eXYwf/wIajfNJ9GrV7mHXrh0MGfIkYP8wcOva\n/vmpWfN+lixZCbh+bQMp/Hfh4EEoV07loYdkkQkhhLgbAS0b5bk9e/xELCNHA+A3djT6uFinNjkt\nW5G2ci0A3hvW4rtwPn8cOX7b40VHf05c3HdkZGRgNCYSGjqMvn37O7Vr0KARO3duJzMzg0OHYtm6\n9TPHc9279+Thh9s4Hn/9dQyPPfY4+/fv4ccf42nRopXT/lq3bseuXV87Cv+uXTto3brdbWMtajLG\nf4dSU+Hnn5GFe4QQopQ5d+433nxzAe++u5xVq97DZnPuvMXGfkv9+g25fPkSNWven+vGPWC/eyyA\nzWZj9+6ddO/em5497XfVy4u/vz9eXl5cunQRi8XCyZM/U79+A9e/uTsgPf47JOvzCyHE35dfD12n\n08D/TmunLVtV4H5MI0ZiGjGyUMds1qwFOp2OSpUq4efnR2rqdQBWrFjCpk0bUFWV+vUb0r//AH7/\n/VyeHwxu+vHHeKpWrUa1atXo3r0X69a9z4svTkOncy6r3br1ZOfO7dStG0yLFq0ct+S9nVvvK6Ao\nCvfdV5OpU18u1PssiBT+OyR35BNCiNLJZvvz1jT2u9TYC/CtY/w31ahRg/Pnf8dsNmMwGBzbT536\nmYceasDXX8dw7dpVRo4cBoDJZOL77w/Srl1Hp+N26dKNSZMmcOnSRR57bACXL18sMNaiHOOXU/13\n6GaPv2VLGd8XQojS5MSJn7BarVy/fp3MzAwqVsx/hryvbzk6duzC6tXvObbt2bOLJUsWYjabOXBg\nP2vXfuT4euGFKfme7q9cuQp+fn6cOnWSxo2buPx93Snp8d+hH3/UUrcuVK5c5m9qKIQQZUq1atV5\n5ZXpXL58kWeeGZvnLPxbTZw4iWXLFhEWNgQ/vwoEBVVlzpy3iYv7jiZNmlKxYiVH227derJy5TKy\ns7Px8vJy2lfXrj34/fdzeR5zzpx/O2b6t2z5MH369HU61a+qKmPHTqBBg7wnRd6JYr8t76FDh5g4\ncSJz5syhW7duAJw6dYpZs2YBEBwczGuvvQbA6tWriYmJQVEUwsPD6dKlC2lpaUyaNIm0tDR8fX15\n5513qFSpUn6HA1x7W96ICC/atzfQr1/Zv5XmnfCU24sWluTDmeTEmafkpCTcljc6+nN+++0s4eHP\nu3zfRa1U35b3woULfPDBB7Ro0SLX9tmzZxMREcHmzZtJT09n7969XLx4kejoaD766CNWrFjB3Llz\nsVqtrFu3jtatW7Np0yZ69+7NqlUFT/5wpTlzsvnnP4v1kEIIIYTLFGvhDwwMZMmSJY7LIQDMZjOX\nL1+mSRP7uEe3bt2IjY0lLi6OTp06YTAYCAgIoEaNGpw5c4bY2Fh69eqVq60QQghxO48++lip7O0X\nhWId4/fx8XHalpKSQoUKFRyPK1eujNFopFKlSgTccgecgIAAjEYjSUlJju2VK1cmMTGxwOP6+/ui\n07n2ovvbnUbxVJKT3CQfziQnzjwhJxaLheTkG/ZL9QqhsO08SX45sdmgcuUKeV5GmO++XBXUX0VG\nRhIZGZlr2/jx4+nUqVM+r7DLb8pBXtsLOz0hJSWzUO0Ky1PG5e6E5CQ3yYczyYkzT8qJ2ZyDxWJF\no7l9J6yoxvhLs/xyYrNZsdlspKRkOT13uw+URVb4Q0JCCAkJKbBdQEAA169fdzxOSEggKCiIoKAg\nzp07l+d2o9GIn5+fY5sQQoiSzdvb11GobqdiRV+Sk9OLKarSIb+c6HT6Aj9I5cXtl/Pp9Xpq1arF\n4cOHadWqFTt27GDEiBE88MADfPDBB4wfP56UlBQSExOpU6cOHTp0ICYmhrFjx7Jjx44CzyAIIYQo\nGTQabYGFymAwFGr2vydxdU6KtfDv2bOHNWvW8Ntvv3HixAk2bNjA+++/T0REBK+++io2m42mTZvS\nvn17AEJDQxk+fDiKojBr1iw0Gg0jRoxgypQpDBs2jAoVKvD2228X51sQQgghSrViv47fHVw9huZJ\n43KFJTnJTfLhTHLiTHLiTHLi7G5yUmKu4xdCCCGEe3lEj18IIYQQdtLjF0IIITyIFH4hhBDCg0jh\nF0IIITyIFH4hhBDCg0jhF0IIITyIFH4hhBDCg0jhF0IIITyI29fqL03mzJnD0aNHURSFiIgImjRp\n4u6Qitzp06cZO3YsI0eOZPjw4Vy9epWpU6ditVoJDAzk7bffxmAwsG3bNtatW4dGoyE0NJSQkBBy\ncnKYPn06V65cQavVMnfuXO677z53v6W/bd68eRw5cgSLxcKzzz5L48aNPTYnWVlZTJ8+neTkZLKz\nsxk7diwPPfSQx+bjViaTiX79+jF27FjatWvn0TmJi4tj4sSJ1K1bF4B69erx9NNPe3ROALZt28bq\n1avR6XRMmDCB4ODg4smJKgolLi5OfeaZZ1RVVdUzZ86ooaGhbo6o6GVkZKjDhw9XZ8yYoW7YsEFV\nVVWdPn26Gh0draqqqr7zzjvqxo0b1YyMDLV3797qjRs31KysLLVv375qSkqK+sknn6izZs1SVVVV\n9+/fr06cONFt78VVYmNj1aefflpVVVX9448/1C5dunh0Tr788kt15cqVqqqq6qVLl9TevXt7dD5u\ntWDBAnXgwIHqxx9/7PE5OXjwoDp+/Phc2zw9J3/88Yfau3dvNS0tTU1ISFBnzJhRbDmRU/2FFBsb\nS8+ePQGoXbs2qamppKeX7VtHGgwGVq1alevWx3FxcfTo0QOAbt26ERsby9GjR2ncuDF+fn54e3vT\nokUL4uPjiY2NpVevXgC0b9+e+Ph4t7wPV3r44Yd59913AahQoQJZWVkenZNHH32U0aNHA3D16lWq\nVq3q0fm46ezZs5w5c4auXbsC8neTF0/PSWxsLO3ataN8+fIEBQXx+uuvF1tOpPAXUlJSEv7+/o7H\nAQEBGI1GN0ZU9HQ6Hd7e3rm2ZWVlYTAYAKhcuTJGo5GkpCQCAgIcbW7m5tbtGo0GRVEwm83F9waK\ngFarxdfXF4CoqCg6d+7s8TkBGDp0KJMnTyYiIkLyAbz11ltMnz7d8VhyAmfOnOG5557jiSee4MCB\nAx6fk0uXLmEymXjuuecYNmwYsbGxxZYTGeO/S6rc4iDfHNzp9tJo586dREVF8f7779O7d2/Hdk/N\nyebNmzl58iRTpkzJ9Z48MR///e9/adasWb7jrZ6YkwceeIDw8HAeeeQRLl68SFhYGFar1fG8J+YE\n4Pr16yxZsoQrV64QFhZWbH870uMvpKCgIJKSkhyPExMTCQwMdGNE7uHr64vJZAIgISGBoKCgPHNz\nc/vNsyI5OTmoqur4NFua7d+/n+XLl7Nq1Sr8/Pw8OifHjx/n6tWrANSvXx+r1Uq5cuU8Nh8Ae/bs\nYdeuXYSGhhIZGcmyZcs8+ncEoGrVqjz66KMoikLNmjWpUqUKqampHp2TypUr07x5c3Q6HTVr1qRc\nuXLF9rcjhb+QOnTowPbt2wE4ceIEQUFBlC9f3s1RFb/27ds78rBjxw46depE06ZNOXbsGDdu3CAj\nI4P4+HhatWpFhw4diImJAWD37t20adPGnaG7RFpaGvPmzWPFihVUqlQJ8OycHD58mPfffx+wD4dl\nZmZ6dD4AFi5cyMcff8zWrVsJCQlh7NixHp+Tbdu2sWbNGgCMRiPJyckMHDjQo3PSsWNHDh48iM1m\nIyUlpVj/duS2vHdg/vz5HD58GEVRmDlzJg899JC7QypSx48f56233uLy5cvodDqqVq3K/PnzmT59\nOtnZ2VSvXp25c+ei1+uJiYlhzZo1KIrC8OHD6d+/P1arlRkzZvD7779jMBh48803ueeee9z9tv6W\nLVu2sHjxYh588EHHtjfffJMZM2Z4ZE5MJhMvv/wyV69exWQyER4eTqNGjZg2bZpH5uOvFi9eTI0a\nNejYsaNH5yQ9PZ3Jkydz48YNcnJyCA8Pp379+h6dE7APkUVFRQEwZswYGjduXCw5kcIvhBBCeBA5\n1S+EEEJ4ECn8QgghhAeRwi+EEEJ4ECn8QgghhAeRwi+EEEJ4EFm5T4gyZt68eRw7dozs7Gx+/vln\nmjdvDsCgQYN4/PHHC7WPlStXUq9ePcda83kZMWIEa9euRavVuiLsXBISEvjtt99o166dy/cthKeT\ny/mEKKMuXbrEsGHD2Ldvn7tDuWPbtm3j7NmzvPDCC+4ORYgyR3r8QniQxYsXc+nSJa5cucK0adMw\nmUzMnz8fg8GAyWRi5syZNGzYkOnTp9OyZUvatWvHmDFj6NixIz/99BMZGRmsWLGCqlWrEhwczIkT\nJ3jvvfe4fv06165d4/z587Rp04ZXXnmF7Oxspk2bxuXLl6lWrRparZYOHToQEhLiiCcjI4NJkyZx\n48YNLBYL3bp1o1+/fixcuBBVValUqRJPPvkk//73vzl//jwZGRn069ePf/3rX3zyySd8/fXXKIpC\nQkICtWrVYs6cOfzxxx9MnjwZsC8wNGTIEAYPHuyulAtR4sgYvxAe5tKlS6xfv55GjRpx/fp1Zs2a\nxfr16wkLC2PFihVO7c+ePcvAgQPZuHEj9evX56uvvnJq8/PPP7No0SKioqL45JNPSE1NZdu2bVgs\nFiIjI3n11Vc5cOCA0+u+++47LBYLH330EZs3b8bX15caNWowYMAA+vfvzz//+U/Wr19PUFAQGzZs\nIDIyki+//JJTp04BcOzYMebPn09UVBRXrlxh3759fPXVV9SqVYsNGzbw4YcfOtY+F0LYSY9fCA/T\ntGlTFEUBoEqVKsybN4/s7GzS0tKoWLGiU3t/f3/q1q0LQPXq1bl+/bpTm5YtW6LVatFqtfj7+5Oa\nmsrJkydp3bo1AIGBjcAhkQAAAlNJREFUgbRs2dLpdS1atGDRokVMnDiRLl26EBISgkaTuz8SFxfH\ntWvX+P777wEwm81cuHDB8fqbt0lu3rw5Z8+epUePHnz00UdMnz6dLl26MGTIkLtNlRBlkvT4hfAw\ner3e8fPUqVMZPXo0GzduzHc8/a+T9/KaFpRXG5vNlquI/7Wgg/0OZZ999hlhYWGcOXOGQYMGOfXQ\nDQYD48aNY8OGDWzYsIHPP//ccStkm83mFFft2rX58ssv6d+/P7GxsYwYMSLP9yWEp5LCL4QHS0pK\nom7dulitVmJiYjCbzS7bd61atfjhhx8ASE5O5siRI05tvv32W/bs2UPLli2ZOnUqvr6+JCcnoygK\nFosFsJ9NuDm8YLPZmDt3ruOsw9GjR8nKykJVVeLj4wkODubzzz/n2LFjtG/fnpkzZ3L16lXHvoQQ\ncqpfCI82evRonnrqKapXr86oUaOYOnUqa9eudcm+Bw4cyJ49exgyZAj33nsvrVq1cjoz8P/t3bGx\nRUAUxvFvxiVTgUhABUJjpEKRCgRC2RVRgBLUIFSCUAGkEkXgRU8B790b7f9XwAn32zOzZ4/v+3q/\n3xqGQZZlKY5jeZ6nKIpU17Vs21ZVVdq2TUVR6DxPpWn6rEQOw1BN02jfdwVBoDiOta6r2raV4zi6\n71tlWer14qgDfjHOB+ArjuPQsizKskzXdSnPc3Vd9/wr8F/jOGqeZ/V9/5F6gCm4BgP4Ctd1NU3T\ns0c8SZKPhT6Av6PjBwDAIDzuAwDAIAQ/AAAGIfgBADAIwQ8AgEEIfgAADPIDH9MvDe9HxLYAAAAA\nSUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 576x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fCB99qb7Xome",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "W_vae, log_sigma_sq = sess.run([vae.dec_weight, vae.log_sigma_sq])\n",
        "W_vae = W_vae.T\n",
        "sigma_sq_vae = np.exp(log_sigma_sq)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "J4DS9pnw8F_m",
        "colab_type": "code",
        "outputId": "5cdc36c5-9d1c-49ea-a3f9-02ef2845011a",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 404
        }
      },
      "source": [
        "avg_log_p_vae = log_p_x_true(W_vae, sigma_sq_vae, DATA_SIZE, covariance)/DATA_SIZE\n",
        "\n",
        "# Reordering rows/columns by eigenvalues for visualization\n",
        "diag = np.diag(W_vae.T.dot(W_vae))\n",
        "sort_idx = np.argsort(-diag)\n",
        "P = np.zeros((W_vae.shape[1], W_vae.shape[1]))\n",
        "for i,j in enumerate(sort_idx):\n",
        "  P[i][j] = 1\n",
        "ordered_WTW = P.dot(W_vae.T.dot(W_vae)).dot(P.T)\n",
        "\n",
        "plt.figure()\n",
        "plt.imshow(ordered_WTW[:30,:30], cmap='gray', vmin=0.0)\n",
        "plt.title(r'Top left block of $W^T W$ (ordered by eigenvalue)')\n",
        "print(\"Avg. VAE Likelihood: {}\".format(avg_log_p_vae))\n",
        "print(\"Avg. Maximum Likelihood: {}\".format(avg_log_p_mle))"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Avg. VAE Likelihood: 941.6113943606014\n",
            "Avg. Maximum Likelihood: 941.6543285914813\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUsAAAFeCAYAAAAMrD22AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAelElEQVR4nO3de1TUdf7H8RdClAQuQoCXzMw0WYHN\nShMvGXgBLFO02lLR0t1Tm7K5Hdd1/SVWdlIp2zK34yVpS+pEy24dKxMys1xXSOkmri3qngpJBYQM\nFYrL9/dHxzkCM/BxBpgZej7O4Rzny/f7mfd7Pvry+52Zz4yPZVmWAAAt6uLuAgDAGxCWAGCAsAQA\nA4QlABggLAHAAGEJAAYISwAw4OfuAtB5bd68Wbm5uSopKVFAQIC6d++u+Ph43Xvvve4uDbhgPrwp\nHe3t//7v/3TnnXfqV7/6lbtLAZzGZTja3ZEjR9S/f393lwG4pFOF5bJly5SYmKjExEQNHjxYcXFx\nttunT592efz8/HyNHz++1f0WLlyoMWPGaNeuXXr99dcveCzT+2nrY885v/6mJk+erJycHNvtL7/8\nUtdcc43eeOMN27ZDhw5pyJAh+v777yVJZ86cUWBgoMvjnO/VV1/Vn/70J+eblLRv3z7Fx8e7NIYz\n47fFHNnzxRdfaO7cuW0+rrNc7XP16tV65pln2rAi13Sq5ywfffRR25/j4+OVnp6uG264ocPreOed\nd5STk6PevXvrxhtv1J133tnhNbjiXP1XXHFFs98FBQU1+o/npZdearZt8+bNmjJlirp166Zjx46p\nR48eLo9zvqNHj2rDhg16++23Xeqzs4mJidGmTZvcXUabSU1N1W233abx48dr8ODB7i6nc51Ztubd\nd9/VrbfeqsTERM2aNUvffPONpJ/+B5w0aZJWrlyphIQExcfH67PPPmtxrO3bt2vSpEkaO3as5syZ\no4qKCklSSkqKGhoaNHfuXA0ePFhVVVVKTExUcXGx3XFWrVqlhIQEJSYm6pNPPjGuWZLefPNNJSQk\nKCEhQX/84x/1448/Njq2trZWKSkpysjIMB73/Po//PDDZscFBgbaAq2iokI7duzQ1KlTVVVVJUn6\n/vvv9dZbb2nmzJmSfjo7vPrqq10e53wvvPCCpk6dajtbddTL0aNHNWrUKD3xxBO2cZ5//nmNGTNG\nU6ZM0b///W/bmI7m094YjvZtaXx7ms79tGnTtG3bNtvvP/jgA02ePNnusfZqaHomt27dOsXGxmra\ntGl65ZVXGp3l2jv+XK8vv/yyJk2apNGjR2vr1q2SpNtvv73RlcD27dttJwF///vflZSUpAkTJmjG\njBkqKSlpVGvTupredvR4+vv7a9asWVq3bl2Lj2OHsTqpuLg4a+/evbbbJSUl1vXXX2999dVXlmVZ\n1qZNm6zZs2dblmVZeXl5VmRkpPXOO+9YlmVZr7/+ujV58uRmY+bl5Vnjxo2zvvnmG2vIkCHWf//7\nX8uyLGvdunVWamqqbb+BAwdax44ds4qLi63IyEi79Z27z7ffftuyLMvKysqy3ee5+2mp5uLiYmv4\n8OHW8ePHrYaGBmvevHnWxo0bbcdalmWlpaVZS5cubXbfLY17fv32LFy40Fq7dq1lWZb117/+1Xr0\n0UetF1980Vq1apVlWZaVkZFhzZkzx7b/pk2brOzsbJfHOV9sbKz1n//8p9VeiouLrcGDB1v//Oc/\nLcuyrEOHDllDhw61ysrKrLq6OuuBBx6w4uLiWpzPpmO0tK+j8ZtyNPcZGRnWvHnzbPv9+c9/ttav\nX9/seEc1nD/3RUVF1vXXX2+dOHHCqqmpsWbOnGmrxdHxxcXF1i9/+Utr8+bNlmVZ1tatW63x48db\nlmVZGzZssBYtWmSrYdGiRVZGRoZVXl5uRUVF2f6+LF682FqyZImtz3HjxjWq6/ztrT2elmVZpaWl\nVlRUlHX27Nlmj0NH+9mcWe7evVs33nij+vbtK0m64447lJ+fr7q6OklSQECAkpKSJEkTJkzQwYMH\nVV1dbXesjz76SMOGDdPAgQMlSXfddZd27Nih+vr6C6rp4osvtt1nUlKSDh48qB9++MGo5t27d2vI\nkCGKiIiQj4+PVq9erXvuucd27KuvvqpvvvlGaWlpF/xYtOTcpXJtba1ee+01paSk6NJLL9Xp06fV\n0NCgV155RbNmzbLtP2fOHE2bNs3lcc45evSoqqqqdM011xj1UltbazuL2bt3r4YOHarLLrtMvr6+\nuu222yS1Pp/nj9HSvo7Gt8fe3E+cOFG7du1SVVWV6uvr9cEHH9j2OZ+jGhoaGmz77N27V8OGDVN4\neLguvvjiRnPQ0vF1dXWaOnWqJGnw4MH69ttvJUmJiYn68MMPVV9fr7q6Ou3cuVOJiYkKDQ1VQUGB\n7amWG264weFVlD2tPfZhYWEKCwvTgQMHjMdsL53qOcuWVFZWNnruKygoSJZlqbKyUpLUrVs3+fj4\n2P4s/XQp2LVr12ZjVVVVad++fUpMTLRtCwwM1HfffafQ0FDjmoKDg9WlSxfb8ZJ06tQpo5qb/u7i\niy+2/bm8vFyrV69WfHy8/PyaT3FL44aFhbVYc1BQkCorK5WTk6NBgwapX79++vLLL1VVVaWdO3fK\n19dXN910U6u9OztORUVFo8ettXn19fVt9NgGBQXZ9j13XEvz2XSMlvZ1NL499ubex8dHMTExys3N\n1RVXXKHevXurT58+zY51VMO5nqWf/u7+4he/sN2OiIho9fjvvvtOvr6+CggIkCR16dLFFsB9+vRR\nz5499emnn6q2tlb9+vVTz549VV9frzVr1tgC7syZM+rXr5/Dvk17Of/fUkhISKOnOtzlZxOWoaGh\n+vTTT223T506pS5duqh79+6SZPuHce530k9/oe0JDw/XiBEjtGbNGpdqOj8Yz73ie/59tlRz9+7d\nG/3u9OnTqqmpkfTTcz1vvPGGZs+erffee6/ZK5KtPRYtCQwMVHFxsV5++WX9/ve/t207ffq0MjMz\nlZKSYvtPpz3GsZq8LbilXo4fP95o327dutmeE5VkC5eW5vPo0aONbre0r6Px7XE097fccou2bdum\nvn37auLEiXaPdVRDfn6+7c+BgYE6e/as7XZpaWmrxzfttamEhAS9//77qq2ttZ3xbt26VTt27FBm\nZqZCQkL0+uuv66233mp0nK+vb6OrrvPf3dBW/5Y6ws/mMnzkyJHat2+f7RLhtdde08iRI21nXjU1\nNdq+fbskKScnR1FRUY3O1s43atSoRmN98cUXevzxx5vtd9FFF6mhocHh25Zqamr03nvv2e4zOjpa\n/v7+RjWPGTNGn3zyiY4ePSrLsrRs2TJlZ2dL+ukfba9evbRixQo9+uijzf5Xbu2xaElQUJD27t2r\nM2fOaNSoUZJ++od58OBB7d+/X8nJya2O4co4ISEh+u6772xnPBfSy5AhQ1RQUKCKigrV19dry5Yt\nkszns7V9HY1vj6O5T0xMVEFBgbZt22b3Ety03piYGOXn56uiokI//vij3nzzzQs63p6EhATt2bNH\nH3zwge1M8OTJk+rdu7dCQkJUWVmpd999V2fOnGl0XFhYmMrKynTy5EnV19c3ClOTWioqKhQSEtJq\nfe3tZ3Nm2aNHDz3++ON64IEHVFtbq8svv1zLly+3/b53794qKCjQk08+qdra2hbf3xUeHq7ly5dr\n3rx5qq2t1aWXXqolS5Y02y8sLEzXX3+94uLitH79el133XWNfn/VVVfp008/1erVq9WlSxetXLnS\nuOYePXroscce0+zZs+Xr66vo6Gjde++9jV7Fv+GGG3TLLbfokUceafQ/d2uPRUuCgoJUWlqqefPm\n2bYFBgaqrKxMs2fP1qWXXtqu41x++eUKDAxUUVGRBg0adEG9REZG6q677lJycrLtLK6oqMh4PqWW\n597R+PY4mvvg4GANHTpUp06dUs+ePS+ohtraWts+MTExSk5OVnJysnr27KmJEyfqb3/7W6s9tKRf\nv35qaGhQRESE7bL+1ltv1TvvvKPx48erT58+WrBggX73u99p5cqViouLkyT17dtX06ZN05QpU9Sr\nVy9NnjxZBw8eNKrl5MmTKisr84i3DrHcUT9dvjz88MO2/+nh2dLS0hQeHq758+e7u5R28cgjj2jA\ngAGaMWOGS+NYlmV7KmPnzp165plnGp1heoOsrCx9+OGHev75591dys/nMhydx29/+1v94x//aHa5\n1xl89dVX+uijj1p8Jd1ERUWFhg8frpKSElmWpXfffVfXXnttG1XZMWpra/XSSy/p/vvvd3cpkghL\neKE+ffroN7/5jdHzbN7k2Wef1Zw5c7R06dJGr6o7IyQkRAsWLNA999yjhIQEnTp1SqmpqW1UacdY\nu3atxo8fr5iYGHeXIonLcAAwwpklABggLAHAAGEJACbcsyTdPkle97N//36310Af9OGpP97WR0uc\nfoHniSee0Oeffy4fHx8tWbKkTV6xMlkm52ms897L5s3ow7PQh3u0FIdOreD5+OOP9fXXXysrK0tH\njhzRkiVLlJWV5XSBAODpnHrOcs+ePRo3bpwkqX///jp16lSbfG0DAHgqp8KyvLy80SfUhISEqKys\nrM2KAgBP0yYfpOHk057tNk5H89a6m6IPz0IfnsWpsAwPD1d5ebntdmlpaasfGmvCm54IPsfbnsB2\nhD48C324R0vB7tRl+MiRI21fXnTgwAGFh4c3+6pTAOhMnDqzvO666zR48GDddddd8vHx0bJly9q6\nLgDwKB71QRredLp+jrddZjhCH56FPtyjzd9neaEcfZdK0+2m79X89a9/3SZ1AYAp1oYDgAHCEgAM\nEJYAYICwBAADhCUAGCAsAcAAYQkABghLADBAWAKAgQ5Z7mi6gsfUK6+8YrTfjBkznBr/Qnjbci5H\n6MOz0Id7tPmnDgHAzw1hCQAGCEsAMEBYAoABwhIADBCWAGCAsAQAA4QlABggLAHAgFeu4DG1cuVK\no/0WL17s9H142woFR+jDs9CHe7CCBwBcRFgCgAHCEgAMEJYAYICwBAADhCUAGCAsAcAAYQkABghL\nADDQqVfwmEpLSzPa77HHHmu2zZP6cAV9eBb6cA9W8ACAiwhLADBAWAKAAcISAAwQlgBggLAEAAOE\nJQAYICwBwABhCQAGWMFzAebPn99s23PPPafU1NRG255//nmj8RoaGtqkrrbgjfNhD314Fm/ro6U4\n9HNmwPz8fD344IMaMGCAJGngwIFaunSpc9UBgBdwKiwladiwYVqzZk1b1gIAHovnLAHAgNNhefjw\nYd1///26++67tXv37rasCQA8jlMv8Jw4cUIFBQVKSkpScXGxZs2apdzcXPn7+9vdv7CwUFFRUS4X\nCwDu0iavht9+++36y1/+oj59+ti/E14Nt4tXw9sefXgWb+ujzT/PcsuWLdq0aZMkqaysTCdPnlRE\nRIRz1QGAF3Dq1fD4+HgtXLhQ77//vmpra/XII484vAQHgM7AqbAMDAzUunXr2roWAPBYrOC5AF26\nNH/Wor6+Xr6+vo223XHHHUbjZWVltUldbcEb58Me+vAs3tYH38EDAC4iLAHAAGEJAAYISwAwQFgC\ngAHCEgAMEJYAYICwBAADhCUAGGAFj4tc6eO2224z2m/Lli1OjX8hmA/PQh/uwQoeAHARYQkABghL\nADBAWAKAAcISAAwQlgBggLAEAAOEJQAYICwBwAAreFzUEX0kJiYa7bdt2zan74P58Cz04R6s4AEA\nFxGWAGCAsAQAA4QlABggLAHAAGEJAAYISwAwQFgCgAHCEgAMsILHRZ7Ux6hRo4z2+9e//tVsmyf1\n4Qr68Cze1gcreADARYQlABggLAHAAGEJAAYISwAwQFgCgAHCEgAMEJYAYICwBAADfu4uAG3H3soc\ne2688Uaj7fn5+S7XBHQWRmeWRUVFGjdunDIzMyVJx44dU0pKiqZPn64HH3xQP/74Y7sWCQDu1mpY\nnj17VsuXL1dsbKxt25o1azR9+nS9+uqr6tu3r7Kzs9u1SABwt1bD0t/fXxs3blR4eLhtW35+vsaO\nHStJiouL0549e9qvQgDwAK0+Z+nn5yc/v8a7VVdXy9/fX5IUGhqqsrKyFsfYv3+/oqKimm3vgA88\n6hCdpY+8vDx3l9AmOst80IdncfkFHpMHIjo62u5x3vTRTY54Yx/2XuDJy8vT8OHDG23zxhd4vHE+\n7KEP92jzj2gLCAhQTU2NJOnEiRONLtEBoDNyKixHjBihnJwcSVJubq5Gjx7dpkUBgKdp9TK8sLBQ\nq1atUklJifz8/JSTk6OnnnpKixcvVlZWlnr16qUpU6Z0RK0A4DathmVUVJQ2b97cbPuLL77YLgUB\ngCfiO3hc1Jn7uPbaa42O/eyzz9qjJKd05vnwRt7WB9/BAwAuIiwBwABhCQAGCEsAMEBYAoABwhIA\nDBCWAGCAsAQAA4QlABhgBY+L6EOKiYkx2u+LL75wavwLwXx4Fm/rgxU8AOAiwhIADBCWAGCAsAQA\nA4QlABggLAHAAGEJAAYISwAwQFgCgAFW8LiIPswNGjTIaL8vv/zS6ftgPjyLt/XBCh4AcBFhCQAG\nCEsAMEBYAoABwhIADBCWAGCAsAQAA4QlABggLAHAACt4XEQfbW/gwIGt7lNUVGR3uyf14Qr6cA9W\n8ACAiwhLADBAWAKAAcISAAwQlgBggLAEAAOEJQAYICwBwICfuwsAmnL0hvPz9e/f3+h3R44caZOa\nAM4sAcCAUVgWFRVp3LhxyszMlCQtXrxYkyZNUkpKilJSUrRz5872rBEA3K7Vy/CzZ89q+fLlio2N\nbbT9oYceUlxcXLsVBgCepNUzS39/f23cuFHh4eEdUQ8AeKRWzyz9/Pzk59d8t8zMTL344osKDQ3V\n0qVLFRIS4nCM/fv3Kyoqqtn2DvjAow5BH57l8OHD7i6hTXSW+egsfTj1avjkyZMVHBysyMhIbdiw\nQWvXrlVaWprD/aOjo5tt87aPbnKEPtzD0avhhw8f1tVXX2277a2vhnvbfDjibX20+Ue0xcbGKjIy\nUpIUHx9v9FYPAPBmToVlamqqiouLJUn5+fkaMGBAmxYFAJ6m1cvwwsJCrVq1SiUlJfLz81NOTo5m\nzpypBQsWqGvXrgoICNCKFSs6olYAcBu+VsJF9OFZmvZx1VVXGR33v//9r71KckpnnQ9Px9dKAICL\nCEsAMEBYAoABwhIADBCWAGCAsAQAA4QlABggLAHAAGEJAAZYweMi+vAszvZx5ZVXGu331VdfXfDY\nzvi5z4e7sIIHAFxEWAKAAcISAAwQlgBggLAEAAOEJQAYICwBwABhCQAGCEsAMMAKHhfRh2dp7z76\n9OljtN+5bz91FvPhHqzgAQAXEZYAYICwBAADhCUAGCAsAcAAYQkABghLADBAWAKAAcISAAywgsdF\n9OFZPKWPyy+/3Gi/o0eP2t3uKX24ytv6YAUPALiIsAQAA4QlABggLAHAAGEJAAYISwAwQFgCgAHC\nEgAMEJYAYMDP3QUAnZGjlTlNtbTS5/zfmY6H9mMUlunp6SooKFBdXZ3uu+8+RUdHa9GiRaqvr1dY\nWJiefPJJ+fv7t3etAOA2rYZlXl6eDh06pKysLFVWVio5OVmxsbGaPn26kpKS9PTTTys7O1vTp0/v\niHoBwC1afc5y6NChevbZZyVJ3bp1U3V1tfLz8zV27FhJUlxcnPbs2dO+VQKAm7Ualr6+vgoICJAk\nZWdn66abblJ1dbXtsjs0NFRlZWXtWyUAuJnxCzzbt29Xdna2MjIyNGHCBNt2k094279/v6Kioppt\n74BPh+sQ9OFZOksfxcXF7i6hTXSW+TAKy127dmndunV64YUXFBQUpICAANXU1OiSSy7RiRMnFB4e\n3uLx0dHRzbZ52+fcOUIfnsXb+nD0anhxcbH69Olju+2tr4Z723y49HmWVVVVSk9P1/r16xUcHCxJ\nGjFihHJyciRJubm5Gj16dBuVCgCeqdUzy61bt6qyslILFiywbVu5cqUefvhhZWVlqVevXpoyZUq7\nFgkA7sbXSriIPjyLt/XBZbhnaSkOCUsX0Ydn6ax99OrVy+i4b7/9tr1Kcoq3zQffwQMALiIsAcAA\nYQkABghLADBAWAKAAcISAAwQlgBggLAEAAOEJQAYYAWPi+jDs/zc++jRo4fRfsePH7/gsZ3hbfPB\nCh4AcBFhCQAGCEsAMEBYAoABwhIADBCWAGCAsAQAA4QlABggLAHAACt4XEQfnoU+zERERBjtd+LE\nCZfux9vmgxU8AOAiwhIADBCWAGCAsAQAA4QlABggLAHAAGEJAAYISwAwQFgCgAFW8LiIPjwLfbSt\n8PBwo/1KS0vtbveUPkyxggcAXERYAoABwhIADBCWAGCAsAQAA4QlABggLAHAAGEJAAYISwAw4Ofu\nAgB4Lkcrc5oKCwsz+l1ZWZnLNbmLUVimp6eroKBAdXV1uu+++7Rjxw4dOHBAwcHBkqS5c+fq5ptv\nbs86AcCtWg3LvLw8HTp0SFlZWaqsrFRycrKGDx+uhx56SHFxcR1RIwC4XathOXToUMXExEiSunXr\npurqatXX17d7YQDgSS7oU4eysrK0b98++fr6qqysTLW1tQoNDdXSpUsVEhLi+E741CGPRx+exdv6\ncPScZWlpaaNPLvL05yxbikPjsNy+fbvWr1+vjIwMFRYWKjg4WJGRkdqwYYOOHz+utLQ0h8cWFhYq\nKirqwisHAE9hGfjoo4+sadOmWZWVlc1+d+jQIWvGjBktHi+p2Y+j7d72Qx+e9UMf7vkJCwuz+2NZ\nVqPb7q7T5HF3pNX3WVZVVSk9PV3r16+3vfqdmpqq4uJiSVJ+fr4GDBjQ2jAA4NVafYFn69atqqys\n1IIFC2zbpk6dqgULFqhr164KCAjQihUr2rVIAHA3vlbCRfThWejDPXiBp40Qlp6PPjxLZ+3jsssu\nMzquvLy8vUpqUUtxyNpwADBAWAKAAcISAAwQlgBggLAEAAOEJQAYICwBwABhCQAGCEsAMMAKHhfR\nh2ehD8/ibB/uWunDCh4AcBFhCQAGCEsAMEBYAoABwhIADBCWAGCAsAQAA4QlABggLAHAACt4XEQf\nnoU+PEt79xESEmK0X0VFhdF+rOABABcRlgBggLAEAAOEJQAYICwBwABhCQAGCEsAMEBYAoABwhIA\nDPi5uwAAcJbpypzu3bu7fF+cWQKAAcISAAwQlgBggLAEAAOEJQAYICwBwABhCQAGCEsAMEBYAoCB\nDvkOHgDwdpxZAoABwhIADBCWAGCAsAQAA4QlABggLAHAgFs+/PeJJ57Q559/Lh8fHy1ZskQxMTHu\nKMMl+fn5evDBBzVgwABJ0sCBA7V06VI3V2WuqKhIDzzwgO655x7NnDlTx44d06JFi1RfX6+wsDA9\n+eST8vf3d3eZrWrax+LFi3XgwAEFBwdLkubOnaubb77ZvUUaSE9PV0FBgerq6nTfffcpOjraK+ej\naR87duzwyvmwp8PD8uOPP9bXX3+trKwsHTlyREuWLFFWVlZHl9Emhg0bpjVr1ri7jAt29uxZLV++\nXLGxsbZta9as0fTp05WUlKSnn35a2dnZmj59uhurbJ29PiTpoYceUlxcnJuqunB5eXk6dOiQsrKy\nVFlZqeTkZMXGxnrdfNjrY/jw4V43H450+GX4nj17NG7cOElS//79derUKZ0+fbqjy/hZ8/f318aN\nGxUeHm7blp+fr7Fjx0qS4uLitGfPHneVZ8xeH95o6NChevbZZyVJ3bp1U3V1tVfOh70+6uvr3VxV\n2+nwsCwvL2/0fRghISEqKyvr6DLaxOHDh3X//ffr7rvv1u7du91djjE/Pz9dcskljbZVV1fbLvNC\nQ0O9Yk7s9SFJmZmZmjVrlv7whz8Yf0eLO/n6+iogIECSlJ2drZtuuskr58NeH76+vl43H464/QvL\nvHW15ZVXXqn58+crKSlJxcXFmjVrlnJzc73ieaXWeOucSNLkyZMVHBysyMhIbdiwQWvXrlVaWpq7\nyzKyfft2ZWdnKyMjQxMmTLBt97b5OL+PwsJCr52Ppjr8zDI8PFzl5eW226WlpQoLC+voMlwWERGh\niRMnysfHR1dccYUuu+wynThxwt1lOS0gIEA1NTWSpBMnTnjtpW1sbKwiIyMlSfHx8SoqKnJzRWZ2\n7dqldevWaePGjQoKCvLa+Wjah7fOhz0dHpYjR45UTk6OJOnAgQMKDw9XYGBgR5fhsi1btmjTpk2S\npLKyMp08eVIRERFursp5I0aMsM1Lbm6uRo8e7eaKnJOamqri4mJJPz0Pe+7dCp6sqqpK6enpWr9+\nve1VY2+cD3t9eON8OOKWTx166qmntG/fPvn4+GjZsmUaNGhQR5fgstOnT2vhwoX6/vvvVVtbq/nz\n52vMmDHuLstIYWGhVq1apZKSEvn5+SkiIkJPPfWUFi9erB9++EG9evXSihUrdNFFF7m71BbZ62Pm\nzJnasGGDunbtqoCAAK1YsUKhoaHuLrVFWVlZeu6559SvXz/btpUrV+rhhx/2qvmw18fUqVOVmZnp\nVfPhCB/RBgAGWMEDAAYISwAwQFgCgAHCEgAMEJYAYICwBAADhCUAGCAsAcDA/wMAiAvQrrQE9gAA\nAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 576x396 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ID1E_hLL8Gww",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}